Method and system for approximate quantum circuit synthesis using quaternion algebra

ABSTRACT

This application concerns methods, apparatus, and systems for performing quantum circuit synthesis and/or for implementing the synthesis results in a quantum computer system. In certain example embodiments: a universal gate set, a target unitary described by a target angle, and target precision is received (input); a corresponding quaternion approximation of the target unitary is determined; and a quantum circuit corresponding to the quaternion approximation is synthesized, the quantum circuit being over a single qubit gate set, the single qubit gate set being realizable by the given universal gate set for the target quantum computer architecture.

CROSS REFERENCE TO RELATED APPLICATIONS

This application is the U.S. National Stage of International Application No. PCT/US2016/053808 entitled “METHOD AND SYSTEM FOR APPROXIMATE QUANTUM CIRCUIT SYNTHESIS USING QUATERNION ALGEBRA,” filed on Sep. 26, 2016, which was published in English under PCT Article 21(2), which in turn claims the benefit of U.S. Provisional Application No. 62/233,293 entitled “METHOD AND SYSTEM FOR APPROXIMATE QUANTUM CIRCUIT SYNTHESIS USING QUATERNION ALGEBRA,” filed on Sep. 25, 2015. The provisional application is incorporated herein by reference in its entirety.

This application is the U.S. National Stage of International Application No. PCT/US2016/053808 entitled “METHOD AND SYSTEM FOR APPROXIMATE QUANTUM CIRCUIT SYNTHESIS USING QUATERNION ALGEBRA,” filed on Sep. 26, 2016, which was published in English under PCT Article 21(2), which is also a continuation-in-part of PCT International Application No. PCT/US2016/025958 filed on Apr. 5, 2016, and entitled “METHOD AND SYSTEM FOR QUANTUM CIRCUIT SYNTHESIS USING QUATERNION ALGEBRA”, which claims the benefit of U.S. Provisional Patent Application No. 62/146,182 filed on Apr. 10, 2015, and entitled “METHOD AND SYSTEM FOR QUANTUM CIRCUIT SYNTHESIS USING QUATERNION ALGEBRA”, both of which are hereby incorporated herein by reference in their entirety.

FIELD

The application pertains to quantum circuit synthesis.

SUMMARY

Disclosed herein are example embodiments of methods, apparatus, and systems for performing quantum circuit synthesis and/or for implementing the synthesis results in a quantum computer system. For instance, the disclosure presents innovations in quantum computing design and synthesis tools for generating circuit descriptions operable on such designs. The innovations can be implemented as part of a method, as part of a computing device, compilation system, or synthesis system configured or programmed to perform any embodiment of the disclosed techniques, or as part of a tangible computer-readable media storing computer-executable instructions for causing a computing device to perform any embodiment of the disclosed synthesis techniques. The innovations can also be performed as part of a method or system for implementing and performing quantum operations in a physical implementation of a quantum computer. The various innovations can be used in combination or separately.

Embodiments of the disclosed technology can be used as part of an overall synthesis process for generating gate sets for a quantum computer having a target architecture from higher-level descriptions. In particular embodiments, the techniques disclosed herein can be used to generate a quaternion approximation for a synthesis process that uses quaternion algebra. Example embodiments of suitable quaternion-based synthesis procedures with which the disclosed technology can be used are described in U.S. Provisional Patent Application No. 62/146,182 filed on Apr. 10, 2015, and entitled “METHOD AND SYSTEM FOR QUANTUM CIRCUIT SYNTHESIS USING QUATERNION ALGEBRA”; PCT International Application No. PCT/US2016/025958 filed on Apr. 5, 2016, and entitled “METHOD AND SYSTEM FOR QUANTUM CIRCUIT SYNTHESIS USING QUATERNION ALGEBRA”; Vadym Kliuchnikov and Jon Yard, “A Framework for Exact Synthesis,” arXiv:1504.04350 (2015); and in Simon Forest, David Gosset, Vadym Kliuchnikov, David McKinnon, “Exact synthesis of single-qubit unitaries over Clifford-cyclotomic gate sets,” Journal of Mathematical Physics 56, 082201, (2015) (also available at arXiv:1501.04944v2). For instance, any of the embodiments disclosed in these references (or aspect or feature (including combinations thereof) disclosed therein) can be used together with any aspect or feature (including combinations thereof) of the innovations presented in this disclosure.

In certain example embodiments, methods for efficiently approximating qubit unitaries over gate sets derived from totally definite quaternion algebras are presented. Example implementations of the technology can achieve ε-approximations using circuits of length O(log(1/ε)), which is asymptotically optimal. Example implementations also achieve the same quality of approximation as previously-known algorithms for Clifford+T, V-basis, and Clifford+π/12, running on average in time polynomial in O(log(1/ε)) (conditional on a number-theoretic conjecture).

In certain example embodiments of the quantum circuit synthesis procedure disclosed herein, a universal gate set, a target unitary described by a target angle, and target precision is received (input); a corresponding quaternion approximation of the target unitary is determined; and a quantum circuit corresponding to the quaternion approximation is synthesized, the quantum circuit being over a single qubit gate set, the single qubit gate set being realizable by the given universal gate set for the target quantum computer architecture. The single qubit gate set can be a Clifford+T basis gate set, a Clifford+e^(iπZ/12) basis gate set, a Clifford+e^(iπZ/16) basis gate set, a V-basis gate set, or any other single qubit base set expressible in terms of totally definite quaternion algebra. The corresponding quaternion approximation can describe the target unitary in quaternion algebra. In some implementations, the system further comprises a quantum circuit controller coupled to the target quantum circuit architecture and configured to implement the quantum circuit in the target quantum circuit architecture. In particular implementations, the synthesis procedure further comprises determining a cost vector for the target unitary, and the act of determining the corresponding quaternion approximation comprises finding one or more corresponding quaternion approximations that satisfy the cost vector. The cost vector can, for example, describe a limit on a size of the single qubit gate set resulting from the synthesis procedure. In some implementations, the act of determining the corresponding quaternion approximation of the target unitary comprises selecting a first algebraic integer for use as a first quaternion element in the quaternion approximation, the first algebraic integer being selected so that a distance threshold to a R_(x), R_(y), or R_(z) rotation is satisfied. The act of determining the corresponding quaternion approximation of the target unitary can further comprise selecting a second algebraic integer for use as a second quaternion element in the quaternion approximation, the second algebraic integer being selected so that, in combination with the first integer, the quaternion approximation produces a quaternion that satisfies a cost limit and the corresponding unitary satisfies the distance threshold. The first algebraic integer and the second algebraic integer can be selected from a ring of integers of the CM-field K. In particular implementations, the determining a corresponding quaternion approximation of the target unitary comprises finding a quaternion from an order of a totally definite quaternion algebra defined over totally real number field F (e.g., the generalized Lipschitz order) that has the following two properties: (1) (U_(q),R_(z)(φ)≤ε; and (2) nrd(q)

_(F)=

₁ ^(L) ¹ . . .

_(M) ^(L) ^(M) , where d is a distance function (e.g., d=√{square root over (1−|Tr(U_(q) ^(d)R_(z)(φ)|/2))},

₁, . . . ,

_(M) are the appropriate prime ideals of

_(F), and L₁, . . . , L_(M)∈

are their respective multiplicities in the decomposition of nrd(q)

_(F).

In some example embodiments of a quantum circuit synthesis procedure disclosed herein, a program describing a desired computation to be performed in a target quantum computer architecture is input; and a gate set adapted for implementation on the target quantum computer architecture is generated based on the program, the gate set including a single qubit circuit for one or more respective unitaries used to perform the desired computation. In certain implementations, the single qubit circuit is generated by identifying a target unitary, and determining a corresponding quaternion approximation of the target unitary using a process in which one or more of the elements in the corresponding quaternion approximation are randomly selected. For example, in some implementations, the act of determining the corresponding quaternion approximation of the target unitary can comprise randomly selecting one or more values for use in the corresponding quaternion approximation from a constrained body of values that guarantees that a norm equation is solvable; and solving the norm equation to determine at least another one of the elements in the corresponding quaternion approximation. For instance, the at least another one of the elements in the corresponding quaternion approximation can be z, where z is an algebraic integer from K, and wherein the norm equation for z is z(z*)=e, where e is an element of the totally real subfield F. In some implementations, the act of determining the corresponding quaternion approximation of the target unitary comprises: randomly sampling points from subsets of a ring of integers; and using the randomly selected sampling points as quaternion elements in the quaternion approximation. The ring of integers can be, for instance, a ring of integers of the CM-field K. In certain implementations, the act of determining the corresponding quaternion approximation of the target comprises: partitioning a convex body of complex numbers into convex subsets, the partitioning being performed such that each subset has a lattice point corresponding to

_(K); randomly selecting one of the subsets; selecting a lattice point corresponding to

_(K) from the randomly selected subset; and using the selected lattice point as an element in the quaternion approximation. In some implementations, the circuit is implemented in the target quantum computer architecture.

In certain example embodiments of a quantum circuit synthesis procedure disclosed herein, a target unitary described by a target angle and target precision is received (input); a corresponding quaternion approximation of the target unitary is determined; and the corresponding quaternion approximation is synthesized using a Closest Vector Problem (CVP) technique to select values for the corresponding quaternion approximation and to produce a single qubit circuit for the target unitary, the single qubit circuit being realizable by a quantum computer architecture. In certain implementations, the determining the corresponding quaternion approximation of the target unitary comprises selecting values from a Hermite-Korkine-Zolotarev, Block-Korkine-Zolotarev, Lenstra-Lenstra-Lovasz, or other size-reduced basis as an element in the corresponding quaternion approximation. Further, in some implementations, the determining the corresponding quaternion approximation of the target unitary comprises using a Lenstra-Lenstra-Lovasz reduction technique or other lattice basis reduction technique.

Any of these example synthesis procedures can be performed by a quantum circuit synthesizer system comprising a processor and memory and/or by a quantum computer synthesis tool adapted for use in a quantum computer design and implementation process and implemented by one or more computing devices. Further, any of these example synthesis procedures can be implemented as computer-executable instructions stored on a computer-readable media, which when executed by a computer cause the computer to perform the synthesis procedure.

BRIEF DESCRIPTION OF THE DRAWINGS

FIG. 1 illustrates a description of an example embodiment of the approximate quaternion approach disclosed herein.

FIG. 2 illustrates a description of an example RANDOM-INTEGER-POINT procedure.

FIG. 3 illustrates a description of an example SUITABLE-Q-NORM procedure.

FIG. 4 illustrates a description of an example IS-EASILY-SOLVABLE procedure.

FIG. 5 illustrates a description of an example FAST-SOLVE-NORM-EQ procedure.

FIG. 6 illustrates a description of an example PASS-HASSE-SOLVABILITY-TEST procedure.

FIG. 7 illustrates a description of an example TOTALLY-POS-GENERATOR procedure.

FIG. 8 illustrates a description of an example TARGET-SIZE procedure.

FIG. 9 illustrates a description of an example UNIT-ADJUST procedure.

FIG. 10 illustrates a description of an example NEAREST-PLANE procedure.

FIG. 11 illustrates a description of an example UNIT-ADJUST procedure.

FIG. 12 illustrates a description of an example RANDOM-INTEGER-POINT procedure.

FIG. 13 is a flow chart showing how the example approximation procedure of FIG. 1 can be implemented as part of a overall synthesis process.

FIG. 14 is a visualization of the projection of the set used in examples of the RANDOM-INTEGER-POINT procedure discussed herein.

FIG. 15 is a visualization of the sampling scheme used in examples of the RANDOM-INTEGER-POINT procedure discussed herein.

FIG. 16 is an example ideal principality graph.

FIGS. 17A-17B together show a table of various quantum circuit creation approaches.

FIG. 18 is a flow chart showing an example synthesis process for producing a gate set in accordance with embodiments of the disclosed technology.

FIG. 19 is a flow chart showing another example synthesis process for producing a gate set in accordance with embodiments of the disclosed technology.

FIG. 20 is a flow chart showing another example synthesis process for producing a gate set in accordance with embodiments of the disclosed technology.

FIG. 21 is a block diagram showing a representative computing environment in which aspects of the disclosed technology can be implemented.

FIG. 22 is a block diagram showing an exemplary system for implementing the disclosed technology, including quantum processing units and quantum controllers for implementing circuit created using embodiments of the disclosed technology.

DETAILED DESCRIPTION 1. General Considerations

As used in this application, the singular forms “a,” “an,” and “the” include the plural forms unless the context clearly dictates otherwise. Additionally, the term “includes” means “comprises.” Further, the term “coupled” does not exclude the presence of intermediate elements between the coupled items. Further, as used herein, the term “and/or” means any one item or combination of any items in the phrase. Still further, as used herein, the term “optimiz*” (including variations such as optimization and optimizing) refers to a choice among options under a given scope of decision, and does not imply that an optimized choice is the “best” or “optimum” choice for an expanded scope of decisions.

The systems, apparatus, and methods described herein should not be construed as limiting in any way. Instead, the present disclosure is directed toward all novel and nonobvious features and aspects of the various disclosed embodiments, alone and in various combinations and sub-combinations with one another. The disclosed systems, methods, and apparatus are not limited to any specific aspect or feature or combinations thereof, nor do the disclosed systems, methods, and apparatus require that any one or more specific advantages be present or problems be solved. Any theories of operation are to facilitate explanation, but the disclosed systems, methods, and apparatus are not limited to such theories of operation.

Although the operations of some of the disclosed methods are described in a particular, sequential order for convenient presentation, it should be understood that this manner of description encompasses rearrangement, unless a particular ordering is required by specific language set forth below. For example, operations described sequentially may in some cases be rearranged or performed concurrently. Moreover, for the sake of simplicity, the attached figures may not show the various ways in which the disclosed systems, methods, and apparatus can be used in conjunction with other systems, methods, and apparatus. Additionally, the description sometimes uses terms like “produce” and “provide” to describe the disclosed methods. These terms are high-level abstractions of the actual operations that are performed. The actual operations that correspond to these terms will vary depending on the particular implementation and are readily discernible by one of ordinary skill in the art.

2. Introduction to the Disclosed Technology

Each time a new computing device is built, one typically asks the question: what problems can it solve? The same question applies to the quantum computers that are now being developed. When addressing such questions, one usually starts with a crude analysis, asking how resources like time, memory, cost, and the size of the computer scale with the problem size. In particular, how do these resources depend on the particular gate set supported by a quantum computer? The algorithm of Solovay and Kitaev shows that any two universal gate sets are equally good from the perspective of polynomially-scaling resources. See, e.g., Alexei Kitaev, Alexander Shen, and Mikhail Vyalyi, “Classical and quantum computation,” volume 47 of Graduate studies in mathematics, American Mathematical Society (2002); Christopher M. Dawson and Michael A. Nielsen, “The Solovay-Kitaev algorithm,” Quantum Information and Computation, 6(1):81-95 (2005). However, once a small quantum computer is built, more refined questions are asked: How large of a problem can one solve on it? How can one compile our algorithms in the most resource-efficient way possible?

Typically, a circuit implementing a quantum algorithm uses a large number of gates, or local unitaries. Each local unitary must be compiled into the gate set supported by a target fault-tolerant quantum computer. Whereas the unitary groups are uncountable, most promising quantum computer architectures known today (topological or based on error correcting codes) natively support only a finite set of unitary gates. The problem of optimal compilation into circuits over such a gate set can be naturally formulated as that of approximation in such groups.

This disclosure concerns the problem of compiling circuits for single-qubit unitaries, and, more specifically, of approximation by finitely-generated subgroups of SU(2). A systematic description of the problem can be given as follows. Let

⊂SU(2) be a finite set of 2×2 unitary matrices, or gates. Given an arbitrary unitary U∈SU(2), one wants to express it in terms of unitaries from

. In most cases U can not be expressed exactly using elements of

and must be therefore approximated. For a selected absolute precision ε, the task is to find a sequence of gates g₁, . . . , g_(N)∈

(usually called a circuit over

) such that ∥U−g_(N) . . . g₁∥≤ε. If one can approximate any unitary over the gate set

(or in other words, if

generates a dense subgroup

of SU(2)), one calls

universal.

Given that each unitary can be so approximated, one can then ask for the shortest, or least costly, such circuit. A volume argument shows that there exist unitaries U requiring circuits of length at least C log(1/ε), where C is a constant that depends on the gate set

. A natural question to ask is whether there is a matching upper bound (e.g., whether one can approximate any unitary using a circuit of length O(log(1/ε))). To answer this question, one can employ non-trivial mathematical ideas. For example, it was recently shown that such approximations exist if the unitaries in

have entries that are algebraic numbers. See, e.g., Jean Bourgain and Alexander Gamburd, “On the spectral gap for finitely-generated subgroups of SU(2),” Inventiones mathematicae, 171(1):83-121 (2008) (hereinafter “Bourgain and Gamurd 2008”). Known gate sets associated with fault-tolerant quantum computing architectures have this property.

Unfortunately, the result is non constructive. Furthermore, there is no obvious way to make it constructive that would realistically work for even a moderately small precision target ε. The result of Bourgain and Gamurd 2008 implies that a brute-force search can yield approximations saturating the lower bound. In practice, however, the precision of approximation achievable with brute force search is limited to 10⁻⁴ or 10⁻⁵. Ideally, one would like to have an algorithm that finds an ε-approximation of a given unitary with a circuit of length O(log(1/ε)) and, furthermore, one would like the algorithm to run in O(poly(log(1/ε)) time.

Recently such algorithms were found for several gate sets such as Clifford+T, the V-basis, Clifford+R_(z)(π/6) and the braiding of Fibonacci anyons. See, e.g., Neil J. Ross and Peter Selinger, “Optimal ancilla-free Clifford+T approximation of z-rotations,” arXiv:1409.4355 (2014); Peter Selinger, “Efficient Clifford+T approximation of single-qubit operators,” arXiv:1212.6253 (December 2012); Alex Bocharov, Yuri Gurevich, and Krysta M. Svore, “Efficient decomposition of single-qubit gates into V basis circuits,” Physical Review A, 88(1):1-13 (July 2013); Alex Bocharov, Martin Roetteler, and Krysta M. Svore, “Efficient synthesis of probabilistic quantum circuits with fallback,” Physical Review A, 91:052317 (2015) (see also arXiv preprint arXiv:1409.3552); Vadym Kliuchnikov, Alex Bocharov, and Krysta M. Svore, “Asymptotically optimal topological quantum compiling,” Physical Review Letters, 112(14) (April 2014). The question of why it is possible to construct such an algorithm for these gate sets and what general properties such gate sets should have has been an outstanding challenge in the field.

In this disclosure, and in certain example embodiments, a general framework is presented that enables efficient approximation algorithms for entire families of gate sets, instead of for specific examples. Example algorithms are developed in the general setting of gate sets derived from totally definite quaternion algebras. Results of applying an example implementation of the algorithm are also presented for a wide range of gate sets, including Clifford+√{square root over (T)}.

The proof that example implementations of the disclosed techniques terminate and run on average in polynomial time relies on a number-theoretic conjecture that generalizes and refines similar conjectures. The mathematics behind conjectures of this type were recently studied for Clifford+T, V-basis and some other gate sets. See, e.g., Peter Sarnak, “Letter to Aaronson and Pollington on the Solvay-Kitaev Theorem and Golden Gates,” available at http://publications.ias.edu/sarnak/paper/2637 (February 2015) (hereinafter “Sarnak 2015”). Results of experiments with example implementations of the disclosed technology provide indirect evidence that some of results in Sarnak 2015 can be true for a wider range of gate sets. This is related to the “Golden Gates” introduced in Sarnak 2015. This is discussed in more detail in the “Conjecture” Section below. Next, the problem of unitary approximation is presented more formally and a high level overview of example embodiments of the approximation framework is presented.

In summary, disclosed herein are embodiments of a framework for approximate synthesis of single qubit unitary transformations over a universal gate set. The example framework is applicable whenever the gate set is related to totally definite quaternion algebras. Embodiments of the disclosed synthesis procedure run in time that is polynomial in log(1/ε), where ε is the approximation parameter and the output factorizations produced have length O(log(1/ε)). An example implementation of the algorithm was implemented in the computer algebra system Magma and was applied to a wide range of gate sets.

2.1. Ancillae Free Approximation.

Formally, the problem of ancillae free approximation for single qubit gate sets can be stated as follows:

Problem 2.1

(unitary approximation problem in two dimensions, UAP). Given

(1) finite universal unitary gate set

⊂SU(2)

(2) cost function c:

→

⁺ (c:

→{1} corresponds to circuit length)

(3) distance function ρ on the set of unitaries

(4) cost bound function cost_(max):

⁺→

⁺

(5) target unitary U from U_(targ)⊂SU(2)

(6) target precision ε

Find g₁, . . . , g_(N) from G such that ρ(g₁· . . . ·g_(N),U)≤ε and Σ_(k=1) ^(N)c(g_(k))≤cost_(max)(ε).

It can be said that an algorithm solves UAP in polynomial time, if it solves Problem 2.1 for arbitrary unitaries U from U_(targ) and its runtime is polynomial in log(1/ε). Here, one also is allowed to spend arbitrary time on precomputation based on (1)-(4) and store an arbitrary amount of results of the precomputation. The set U_(targ) can be equal to SU(2) or some its subset. For example, it can be the set of all unitaries e^(iφz) for Z being Pauli Z matrix and φ being arbitrary real number.

The hardness of solving UAP and the existence of the solution to it depends on the choice of cost bound function cost_(max). A summary of known algorithms for solving UAP is presented in the table shown in FIGS. 17A and 17B. In particular, the table shows known algorithms for solving the unitary approximation problem (UAP, Problem 2.1). The distance function used is ∥U−V∥, where ∥U∥=½√{square root over (Tr(UU^(†)))}. Cost function c_(canonical) is defined below in the “Exact synthesis results for totally definite quaternion algebras” Subsection (Subsection 3.3) and also discussed in more detail in the “End to end examples of using the framework” Section (Section 6).

In practice, for target precisions 10⁻¹⁰ to 10⁻³⁰, the overhead from using the Solovay Kitaev algorithm can be between one to three orders of magnitude. On the other hand, the methods based on brute force search find the best possible solution, but are frequently limited to precisions 10⁻⁵ or even less because their runtime and required memory scale exponentially with log(1/ε). The methods described in Alex Bocharov, Yuri Gurevich, and Krysta M. Svore, “Efficient decomposition of single-qubit gates into V basis circuits,” Physical Review A, 88(1):1-13 (July 2013); Neil J. Ross and Peter Selinger, “Optimal ancilla-free Clifford+T approximation of z-rotations,” arXiv:1409.4355 (2014); and Peter Selinger, “Efficient Clifford+T approximation of single-qubit operators,” arXiv:1212.6253 (December 2012), together with the methods that are described herein (see the Table of FIGS. 17A-17B at entries lines 3-8) are based on an exact synthesis algorithm and produce results that have a cost bound similar to brute force search but have polynomial runtime subject to some number theoretic conjectures (similar to Conjecture 4.8). In the next subsection, this is discussed in more detail.

2.2. Approximation Methods Based on Exact Synthesis.

FIG. 13 is a flow diagram 1300 showing an example approach that solves the UAP (Problem 2.1) and that is based on the exact synthesis algorithm (see the Table of FIGS. 17A and 17B at lines 3-8). The flow diagram is presented at a high level and should not be construed as limiting, as any one or more of the items/acts shown in the diagram 1300 of FIG. 13 can be omitted, replaced, and/or supplemented with additional items/acts.

One of the focuses of this disclosure is to present example algorithms that work for gate sets described by an arbitrary totally definite quaternion algebra. In this subsection, a discussion is provided regarding what it means for the gate set to be described by totally definite quaternion algebra. The flow diagram 1300 of FIG. 13 is also described in more detail.

The “Basic results and definitions” Section (Section 3) below provides definitions and a more detailed discussion of the mathematical objects discussed herein. The aim of this part is to explain connections between them and to the algorithm presented in this disclosure at a high level. It can be said that the gate set

is described by quaternion algebra if the following list of objects can be specified and related to the gate set.

Definition 2.2.

A quaternion gate set specification is a tuple

F, σ, a, b,

, S

where:

-   -   F is a totally real number field and σ is an embedding of F into     -   a, b are elements of F that define the quaternion algebra

$\left( \frac{a,b}{F} \right)$

-   -    over F     -   is a maximal order of

$\left( \frac{a,b}{F} \right)$

-   -   S={         ₁, . . . ,         _(M)} is a set of prime ideals of F

Using the embedding σ any quaternion q from the quaternion algebra can be mapped to a special unitary U_(q)∈SU(2). The precise construction of this map is discussed in the “Using quaternions to represent unitaries” subsection (Subsection 3.1) below. This map has the following relevant properties: U _(q) ₁ _(q) ₂ =U _(q) ₁ U _(q) ₂ ,U _(q) ^(†) =U _(q*) where q* is the conjugate of q. Also defined is the following closed under multiplication set

_(S) ={q∈

:nrd(q)

_(F)=

₁ ^(L) ¹ · . . . ·

_(M) ^(L) ^(M) ,L _(k) ∈

,L _(k)≥0}, and call (L₁, . . . , L_(M)) the cost vector of q. The meaning of a cost vector is discussed in more detail further in this section and also in the “Exact synthesis results for totally definite quaternion algebras” subsection below. Above nrd(q) is the reduced norm of quaternion and

_(F) is a ring of integers of number field F. The set

_(S) is closed under multiplication because

is closed under multiplication and nrd(q₁q₂)=nrd(q₁)·nrd(q₂).

A simplified set of conditions that desirably holds for the gate set to be described by the quaternion gate set specification is:

(1) There must exist subset

_(Q) of

_(S) such that

={U_(q):q∈

_(Q)}.

(2) The group generated by

must be equal to group {U_(q):q∈

_(S)}.

Condition (1) implies that the group generated by elements from G is a subgroup of {U_(q):q∈

_(S)}. Condition (2) can be checked for a given set

_(Q) using the framework described in Vadym Kliuchnikov and Jon Yard, “A framework for exact synthesis,” arXiv:1504.04350 (2015); U.S. Provisional Patent Application No. 62/146,182 filed on Apr. 10, 2015, and entitled “METHOD AND SYSTEM FOR QUANTUM CIRCUIT SYNTHESIS USING QUATERNION ALGEBRA”; and PCT International Application No. PCT/US2016/025958 filed on Apr. 5, 2016, and entitled “METHOD AND SYSTEM FOR QUANTUM CIRCUIT SYNTHESIS USING QUATERNION ALGEBRA”. A brief overview of results from those works is in the “Exact synthesis results for totally definite quaternion algebras” subsection (Subsection 3.3) below.

One way of checking the condition (2) is to first compute a finite set of quaternions

_(,S) such that every element of

_(S) can be written as a product of elements of

_(,S) and a scalar (using, for example, algorithms from Vadym Kliuchnikov and Jon Yard, “A framework for exact synthesis,” arXiv:1504.04350 (2015); U.S. Provisional Patent Application No. 62/146,182 filed on Apr. 10, 2015, and entitled “METHOD AND SYSTEM FOR QUANTUM CIRCUIT SYNTHESIS USING QUATERNION ALGEBRA”; and PCT International Application No. PCT/US2016/025958 filed on Apr. 5, 2016, and entitled “METHOD AND SYSTEM FOR QUANTUM CIRCUIT SYNTHESIS USING QUATERNION ALGEBRA”). Here,

_(,S) is a set of canonical generators of

_(S). Second, for each q from

_(Q)*, find a representation of U_(q) as a product of elements of

. For all q from

_(,S) one then can define: Circuit(q)=(U ₁ , . . . ,U _(n)), where U _(q) =U ₁ . . .U _(n) ,U _(k)∈

.  (1) One natural way of defining the cost of elements of

_(,S) is

${c(q)} = {\sum\limits_{k = 1}^{n}{{c\left( U_{k} \right)}.}}$ For the other cost function definitions related to the cost vector of the quaternion, see the “Exact synthesis results for totally definite quaternion algebras” subsection below.

To summarize, the following definition is given:

Definition 2.3.

a gate set

is described by quaternion algebra if the following data is defined:

(1) A quaternion gate set specification

F, σ, a, b,

, S

,

(2) A set set of canonical generators

_(,S) of

_(S),

(3) A map Circuit as described by equation (1).

In the “End to end examples of using the framework” Section (Section 6) below, examples of the gate sets described by quaternion algebras (including Clifford+T and V basis) are given and items (1)-(3) are explicitly specified for each example. With this background, an example flow of an embodiment of the disclosed technology can be given in more detail.

Consider items 1302 and 1303 in FIG. 13. For simplicity, this discussion focuses on the case when set S consists of only one prime ideal. The input to item 1302 is L₁,φ,ε. The output of item 1302 is a quaternion q from

_(S) such that its cost is equal L₁ and U_(q) is within distance ε from R_(z)(φ). The fact that the target gate set

is described by quaternion algebra (Definition 2.3) immediately implies that U_(q) can be expressed as a circuit over

. Items 1303 and 1304 construct such a circuit. At item 1303, q is expressed as a product q₁ . . . q_(n) of elements of

_(,S)* and a scalar from F using an embodiment of the exact synthesis algorithm described in Vadym Kliuchnikov and Jon Yard, “A framework for exact synthesis,” arXiv:1504.04350 (2015); U.S. Provisional Patent Application No. 62/146,182 filed on Apr. 10, 2015, and entitled “METHOD AND SYSTEM FOR QUANTUM CIRCUIT SYNTHESIS USING QUATERNION ALGEBRA”; and PCT International Application No. PCT/US2016/025958 filed on Apr. 5, 2016, and entitled “METHOD AND SYSTEM FOR QUANTUM CIRCUIT SYNTHESIS USING QUATERNION ALGEBRA”. At item 1304, a circuit for U_(q) over gate set

is found as a concatenation of circuits for each q_(k).

The following paragraphs discuss item 1301 in more detail. To give some intuition, one can start with the Clifford+T gate set example (analyzed in more details in the “End to end examples of using the framework” Section (Section 6) below. In this case, set S contains precisely one prime ideal

₁ and L₁ is greater or equal to the T-count of the resulting circuit. To ensure that the approximation step succeeds, the input to the algorithm desirably satisfies inequality: L ₁ log(N(

₁))≥4 log(1/ε)+C _(min), where N(

₁) is the norm of

₁. This reproduces the result from Peter Selinger, “Efficient Clifford+T approximation of single-qubit operators,” arXiv:1212.6253 (December 2012) that the T-count scales as 4 log₂(1/ε)+C_(min) because N(

₁)=2 holds for the Clifford+T case. The bound also saturates the lower bound proved in that paper up to an additive constant. In this simple case, at item 1301 of FIG. 13, one just assigns L₁=┌(4 log(1/ε)+C_(min))/log(N(

₁))┐. In certain example implementations of the disclosed technology, constant C_(min) is precomputed based on a quaternion gate set specification.

More generally, the cost vector (L₁, . . . , L_(M)) that is input to item 1302 desirably satisfies the following inequality: L ₁ log(N(

₁))+ . . . +L _(M) log(N(

_(M)))≥4 log(1/ε)+C _(min). The length of the circuit output by the example embodiment is proportional to L₁+ . . . +L_(M) and therefore proportional to log(1/ε) which is up to multiplicative factor is the best possible. Cost optimality up to an additive constant is more subtle and is dependent on the choice of cost function and the gate set.

In some cases, one might not have very fine control of the cost of the output circuit using cost vector. In the worst case, the cost vector will allow one to control the cost of the output circuit up to multiplicative factor. In this situation, one can use the following strategy to improve the cost of output

-   -   In Step 1 (item 1301 in FIG. 13) make a request for several cost         vectors     -   In Step 2 (item 1302 in FIG. 13), for each cost vector, request         several quaternions that achieve required quality of         approximation     -   In Step 4 (item 1304 in FIG. 13) try to optimize final circuits         using rewriting rules for gate set         and pick the ones with minimal cost. Rewriting rules for the         gate set         can be obtained using an embodiment of the exact synthesis         framework (see Vadym Kliuchnikov and Jon Yard, “A framework for         exact synthesis,” arXiv:1504.04350 (2015)); U.S. Provisional         Patent Application No. 62/146,182 filed on Apr. 10, 2015, and         entitled “METHOD AND SYSTEM FOR QUANTUM CIRCUIT SYNTHESIS USING         QUATERNION ALGEBRA”; and PCT International Application No.         PCT/US2016/025958 filed on Apr. 5, 2016, and entitled “METHOD         AND SYSTEM FOR QUANTUM CIRCUIT SYNTHESIS USING QUATERNION         ALGEBRA”.

More formally, the problem being solved in item 1302 of FIG. 13 is the following:

Problem 2.4.

(Quaternion approximation problem, QAP). Given

(1) A quaternion gate set specification

F, σ, a, b,

, S={

₁, . . . ,

_(M)}

,

(2) target angle φ

(3) target precision ε

(4) target cost vector (L₁, . . . , L_(M)) satisfying L ₁ log N(

₁)+ . . . +L _(M) log N(

_(M))−4 log(1/ε)∈[C _(min) ,C _(max)],

where constants C_(min),C_(max) depend only on the quaternion gate set specification.

Find q from the generalized Lipschitz order (see the “generalized Lipschitz order” section) in quaternion algebra

$\left( \frac{a,b}{F} \right)$ such that nrd(q)

_(F)=

₁ ^(L) ¹ . . .

_(M) ^(L) ^(M) and ∥R_(z)(φ)−U_(q)∥≤ε, where ∥U∥=½√{square root over (Tr(UU^(†)))}.

The polynomial time algorithm (in log(1/ε)) for QAP provides a polynomial time algorithm for solving the unitary approximation problem for gate sets that can be described by totally quaternion algebra. The circuit for U_(q) can be found in time polynomial in L₁, . . . , L_(M) using an embodiment of the exact synthesis algorithm from Vadym Kliuchnikov and Jon Yard, “A framework for exact synthesis,” arXiv:1504.04350 (2015); U.S. Provisional Patent Application No. 62/146,182 filed on Apr. 10, 2015, and entitled “METHOD AND SYSTEM FOR QUANTUM CIRCUIT SYNTHESIS USING QUATERNION ALGEBRA”; or PCT International Application No. PCT/US2016/025958 filed on Apr. 5, 2016, and entitled “METHOD AND SYSTEM FOR QUANTUM CIRCUIT SYNTHESIS USING QUATERNION ALGEBRA”. The cost of the resulting circuit is a linear function in L₁, . . . , L_(M) and therefore one can solve UAP with cost bound function that is in 8(log(1/ε)). Next, the basic definitions used to describe the details of the example embodiments are provided.

3. Basic Results and Definitions

3.1. Using Quaternions to Represent Unitaries.

Let F be a totally real number field of degree d. Let σ₁, . . . , σ_(d) be embeddings of F into

. Let

_(F) be a ring of integers of F. Let a, b be two totally negative elements of

_(F). In other words for all k=1, . . . , d one has σ_(k)(a)<0, σ_(k)(b)<0. Now consider a quaternion algebra

$Q = \left( \frac{a,b}{F} \right)$ given by Q={a ₀ +a ₁ i+a ₂ j+a ₃ k:a ₀ ,a ₁ ,a ₂ ,a ₃ ∈F}, where i²=a, j²=b and k=ij=−ji. The fact that a, b are totally negative implies that Q is totally definite quaternion algebra. The conjugate of a quaternion q=a₀+a₁i+a₂j+a₃k is defined as q*=a₀−a₁i−a₂j−a₃k. The reduced norm nrd and reduced trace trd are defined as nrd(q)=qq*,trd(q)=q+q* Let σ=σ₁ be a fixed embedding that one can use to construct unitaries out of quaternions. A homomorphism from quaternion algebra into the algebra of complex 2×2 matrices can be defined as follows:

$\begin{matrix} \begin{matrix} {{h(q)} = {h\left( {a_{0} + {a_{1}i} + {a_{2}j} + {a_{3}k}} \right)}} \\ {{{\sigma\left( a_{0} \right)}I} + {{\sigma\left( a_{1} \right)}\sqrt{\sigma(a)}Z} + {{\sigma\left( a_{2} \right)}\sqrt{\sigma(b)}Y} + {{\sigma\left( a_{3} \right)}\sqrt{{- {\sigma(b)}}{\sigma(a)}}X}} \end{matrix} & (2) \end{matrix}$ Here I, X, Y, Z are the four Pauli matrices. Note that h has additional nice properties: det(h(q))=σ(nrd(q)),Tr(h(q))=σ(trd(q)).

To construct special unitaries out of quaternions, one can use the following mapping:

$\begin{matrix} {{U(q)} = {\frac{h(q)}{\sqrt{\det\left( {h(q)} \right)}} = \frac{h(q)}{\sqrt{\sigma\left( {{nrd}(q)} \right)}}}} & (3) \end{matrix}$ Note that for any non-zero a from F one has the following

${U\left( {\alpha\; q} \right)} = {{\frac{\sigma(\alpha)}{{\sigma(\alpha)}}{U(q)}} = {\pm {{U(q)}.}}}$

Now consider the structure of image of U(q) in more detail. This structure can be expressed using a number field that can be embedded into Q. Let K=F(√{square root over (a)}) be a totally imaginary extension of F of degree 2. Such number fields K are called CM fields. See, e.g., Lawrence C. Washington, “Introduction to Cyclotomic Fields,” (Springer, 1982). This is ensured by the condition that σ_(k)(a)<0, k=1, . . . , d. Let β be an element of K such that β²−a=0. The degree of the field K is 2d and there are 2d embeddings of K into

. Each element of K can be represented as a₀+a₁β where a₀ and a₁ are from F. One can define 2d embeddings of K into

as following: σ_(k,+)(a ₀ +a ₁β)=σ_(k)(a ₀)+iσ _(k)(a ₁)√{square root over (|σ_(k)(a)|)} σ_(k,−)(a ₀ +a ₁β)=σ_(k)(a ₀)−iσ _(k)(a ₁)√{square root over (|σ_(k)(a)|)} Further, the notation σ for σ_(1,+) is used, which is in agreement with σ=σ₁ for elements of F because for a₀ from F one has σ_(k,±)(a₀)=σ_(k)(a₀).

Each element of the quaternion algebra q=a ₀ +a ₁ i+a ₂ j+a ₃ k=(a ₀ +a ₁ i)+(a ₂ +a ₃ i)j can be mapped to two elements of K in the following way: e ₁(q)=a ₀ +βa ₁ ,e ₂(q)=a ₂ +a ₃.  (4) Conversely, the map e₁ ⁻¹ describes an embedding of K into quaternion algebra Q. Note that now homomorphism h(q) can be written as:

$\begin{matrix} {{h(q)} = {\begin{pmatrix} {\sigma\left( {e_{1}(q)} \right)} & {{\sigma\left( {e_{2}(q)} \right)}\sqrt{{\sigma(b)}}} \\ {{- {\sigma\left( {e_{2}(q)} \right)}^{*}}\sqrt{{\sigma(b)}}} & {\sigma\left( {e_{1}(q)} \right)}^{*} \end{pmatrix}.}} & (5) \end{matrix}$ Using this notation one also has: σ(nrd(q))=|σ(e ₁(q))|²+|σ(b)∥σ(e ₂(q))|². Or in other words, in terms of relative norm N_(K/F) one has: nrd(q)=N _(K/F)(e ₁(q))−bN _(K/F)(e ₂(q)). For any CM field, one can define an automorphism *:K→K which is called complex conjugation and which has the following properties: σ_(k,±)((a ₀ +a ₁β)*)=σ_(k,±)(a ₀ +a ₁β)*, (a ₀ +a ₁β)*=a ₀ −a ₁β. Using it one can express the relative norm N_(K/F)(x)=xx* and see that σ_(k)(N_(K/F)(x))=|σ_(k,±)(a₀+a₁β)². In addition one has that e ₁(q)*=e ₁(q*),e ₁ ⁻¹(q)*=e ₁ ⁻¹(q*). 3.2. Distance to R_(z) Rotations. Example distance functions that are used for unitaries in example embodiments include

$\begin{matrix} {{{\rho\left( {U,V} \right)} = \sqrt{1 - \frac{{{Tr}\left( {UV}^{\dagger} \right)}}{2}}},{{d_{2}\left( {U,V} \right)} = {\frac{1}{2}{\sqrt{{Tr}\left( {\left( {U - V} \right)\left( {U - V} \right)^{\dagger}} \right)}.}}}} & (6) \end{matrix}$ The notation for R_(z) is the following:

$\begin{matrix} {{R_{z}(\varphi)} = {e^{{- \frac{i\;\varphi}{2}}Z} = {\begin{pmatrix} e^{{- i}\;{\varphi/2}} & 0 \\ 0 & e^{i\;{\varphi/2}} \end{pmatrix}.}}} & (7) \end{matrix}$

Now consider the distance between R_(z)(φ) and unitary U_(q) for a given quaternion q:

${\rho\left( {{R_{z}(\varphi)},U_{q}} \right)} = {{\rho\left( {{R_{z}(\varphi)},\frac{h(q)}{R}} \right)}.}$ where R=√{square root over (σ(nrd(q)))} and h(q) is defined by Equation 2. One can further rewrite this as:

${\rho\mspace{11mu}\left( {{R_{z}(\varphi)},{\frac{1}{R}\begin{pmatrix} z & {- w^{*}} \\ w & z^{*} \end{pmatrix}}} \right)},$ where z=σ(e₁(q)), w=−σ(e₁(q))*√{square root over (|σ(b)|)} and |x|²+|y|²=R². Now, one can solve the inequalities:

$\begin{matrix} {{\rho\mspace{11mu}\left( {{R_{z}(\varphi)},{\frac{1}{R}\begin{pmatrix} z & {- w^{*}} \\ w & z^{*} \end{pmatrix}}} \right)} \leq ɛ} & (8) \\ {{d_{2}\left( {{R_{z}(\varphi)},{\frac{1}{R}\begin{pmatrix} z & {- w^{*}} \\ w & z^{*} \end{pmatrix}}} \right)} \leq ɛ} & (9) \end{matrix}$ Inequalities above do not constrain w. Introducing z₀=R(1−ε²)e^(−iφ/2) inequality (8) simplifies to the following two inequalities: Re((z−z ₀)e ^(iφ/2))≥0 or Re((z+z ₀)e ^(iφ/2))≤0, and inequality (9) simplifies to Re((z−z ₀)e ^(iφ/2))≥0. In addition, the fact that U_(q) is a unitary matrix implies that |z|≤R.

FIG. 14 is a visualization 1400 of the solution regions. The result can be summarized in the following proposition:

Proposition 3.1.

Let q=e₁ ⁻¹ (z₁)+e₂ ⁻¹(z₂) and let σ_(1,+)(z₁) belong to the following set {z∈

:Re((z−z ₀)e ^(iφ/2))≥0,|z|≤R} (region 1412 on FIG. 14) where z₀=R(1−ε²)e^(−iφ/2) and R=√{square root over (σ₁(nrd(q)))}. Then d₂(U_(q),R_(z)(φ))≤ε and ρ(U_(q),R_(z)(φ))≤ε.

More specifically, in FIG. 14, the colored regions 1410 (red region) and 1412 (blue region) correspond to complex numbers z that satisfy inequality (8) and |z|≤R given parameters φ,ε,R. The blue region 1412 corresponds to complex numbers z that satisfy inequality (9) and |z|≤R. The vertical axis is Im(z) and the horizontal axis is Re(z).

3.3. Exact Synthesis Results for Totally Definite Quaternion Algebras.

Recall several definitions about rings, ideals, orders, and quaternions that are useful in studying a special case presented in Vadym Kliuchnikov and Jon Yard, “A framework for exact synthesis,” arXiv:1504.04350 (2015), U.S. Provisional Patent Application No. 62/146,182 filed on Apr. 10, 2015, and entitled “METHOD AND SYSTEM FOR QUANTUM CIRCUIT SYNTHESIS USING QUATERNION ALGEBRA”; or PCT International Application No. PCT/US2016/025958 filed on Apr. 5, 2016, and entitled “METHOD AND SYSTEM FOR QUANTUM CIRCUIT SYNTHESIS USING QUATERNION ALGEBRA”, namely the case of totally definite quaternion algebras. Here, Section 2 of Vadym Kliuchnikov and Jon Yard, “A framework for exact synthesis,” arXiv:1504.04350 (2015) is followed. For further references to the literature and other facts used in the following, refer to the reference section of Vadym Kliuchnikov and Jon Yard, “A framework for exact synthesis,” arXiv:1504.04350 (2015).

Let Q be a quaternion algebra over number field F (defined in Section 3.1). A

_(F)-lattice I is finitely generated

_(F)-submodule of Q such that F I=Q. In other words I has a full rank in Q. An order 0 is a

_(F)-lattice that is a subring of Q. An order is a maximal order if it is not properly contained in any other order. There is a right and left order associated with any

_(F) lattice I defined as

_(R)(I)={q∈Q:Iq⊂I},

_(L)(I)={q∈Q:qI⊂I}. When one wants to emphasize that I has particular right and left order, one can call I right-

_(R)(I) fractional ideal or left-

_(L)(I) fractional ideal. A fractional right-

ideal is a normal ideal if the order

is a maximal order. Note that, order

_(R)(I) is maximal if and only if

_(L)(I) is maximal. All normal ideals are invertible. A normal ideal I is principal if I=q

_(R)I for some q from Q.

A normal ideal I is two sided

-ideal if

_(R)(I)=

_(L)(I)=

. The principal two sided

-ideals form a subgroup of the group of all two sided

-ideals (under multiplication). The quotient of the group of all two sided

ideals modulo principal two sided O-ideals is the two-sided ideal class group of

. It is known that the two sided ideal class group of

is always finite. The two sided class number of Q is the size of the two-sided ideal class group of any maximal order of Q. It known that the size of two-sided ideal class group is independent on the choice of maximal order

. Here, the special case of results from Vadym Kliuchnikov and Jon Yard, “A framework for exact synthesis,” arXiv:1504.04350 (2015) are considered for totally definite quaternion algebras with two sided class number 1. The examples considered in this disclosure have this property.

Let disc(q₁, . . . , q₄)=det(trd(q_(i)q_(j))_(i,j=1, . . . , 4)). The discriminant of an order

is

_(F) ideal generated by the set {disc(q ₁ , . . . ,q ₄):q ₁ , . . . ,q ₄∈

}. It turns out that the discriminant always is a square. Its square root is the reduced discriminant denoted by disc(

). All maximal orders in quaternion algebra q has the same discriminant. The reduced norm of

_(F) lattice I is ZF ideal nrd(I) generated by {nrd(q):q∈I}.

The unit group

^(x) of

is {q∈

:nrd q∈

_(F) ^(x)} where

_(F) ^(x) is a unit group of

_(F). For orders in totally definite quaternion algebras, the quotient group

^(x)/

_(F) ^(x) is always finite. The normalizer of order

is the set Normalizer(

)={q∈

:q

q ⁻¹=

}. which is a monoid under multiplication. For totally definite quaternion algebras, the quotient Normalizer(

)/

_(F) (considered as a quotient of two monoids) is finite similarly to

^(x)/

_(F).

One can say that nrd(q) is supported on the set S of primes ideals of

_(F) if:

$\begin{matrix} {{{{nrd}(q)}{\mathbb{Z}}_{F}} = {\prod\limits_{{\mathfrak{p}} \in S}\;{\mathfrak{p}}^{\upsilon{({q,{\mathfrak{p}}})}}}} & (10) \end{matrix}$ Also recall that map T₂: F→

⁺ is defined as:

$\begin{matrix} {{T_{2}(x)} = {\sum\limits_{k = 1}^{d}{{\sigma_{k}^{2}(x)}.}}} & (11) \end{matrix}$

Now all definitions are in place to state the special case of one of the main results of Vadym Kliuchnikov and Jon Yard, “A framework for exact synthesis,” arXiv:1504.04350 (2015) (Theorem 3.18) for totally definite quaternion algebras with two sided class number 1.

Theorem 3.2.

Let Q be a totally definite quaternion algebra over totally real number field F with two sided class number one, let

be a maximal order in Q, let S=(

₁, . . . ,

_(M)) be a finite set of prime ideals of

_(F). There exists set gen_(S)(

) such that every quaternion q from the set

_(S) ={q∈

:nrd(q) is supported on S} can be written as the product q₁ . . . q_(n)q_(rem) where q₁, . . . , q_(n) are from gen_(S)(

) and q_(rem) is from Normalizer(

). If all ideals from S do not divide disc(

) then q_(rem) is from

^(x).

There exist algorithms for deciding if the set gen_(S)(

) is finite and computing it if this is the case. There is also an algorithm for finding factorization q₁ . . . q_(n)q_(r) in time polynomial in log T₂(nrd(q)).

To find the factorization, one can do trial division of q by elements of gen_(S)(

) and greedily reduce values v(

,q) in equation (10) on each step. A trial-division step can be performed until one is left with an element of Normalizer(

). The map U_(q) discussed in the “Using quaternions to represent unitaries” Section (Section 3.1) depends only on q_(r)/F^(x) up to a sign, therefore there is only finitely many possible unitaries U_(q) _(r) . The canonical gate set corresponding to

_(S) is

_(,S) *=gen _(S)(

)∪Normalizer(

)/

_(F)∪{−1}. The main difficulty of the exact synthesis of quaternions and unitaries is computing gen_(S)(

) such that described simple trial division algorithm works. Specific examples illustrating the above definitions are presented below in the “End to end examples of using the framework” Section below.

The canonical cost function for U_(q) can be defined using v(q,

) (see Equation 10) as:

${{c_{canonical}\left( U_{q} \right)} = {\min\limits_{{q^{\prime} \in {\mathcal{M}_{S}:U_{q}}} = U_{q^{\prime}}}{{cost}\left( q^{\prime} \right)}}},{where}$ ${{cost}(q)} = {\sum\limits_{{{\mathfrak{p}} \in S},{{\mathfrak{p}} \nmid {{disc}{(\mathcal{M})}}}}{v\left( {q,{\mathfrak{p}}} \right)}}$ As will be discussed in more detail in the “End to end examples of using the framework” Section (Section 6), the canonical cost function corresponds to the T-count for Clifford+T case and to the V-count for V-basis case. The cost vector (L₁, . . . , L_(M)) of quaternion q is equal to (v(q,

), . . . , v(q,

_(M))). Given the cost vector it is always possible to upper bound C_(canonical)(U_(q)) as:

${{c_{canonical}\left( U_{q} \right)} \leq {{cost}(q)}} = {\sum\limits_{{k = 1},\ldots\mspace{11mu},{M:{{\mathfrak{p}}_{k} \nmid {{disc}{(\mathcal{M})}}}}}{L_{k}.}}$ For the decomposition q=q₁ . . . q_(n)q_(r) described in the Theorem 3.2 one has cost(q)=Σ_(k=1) ^(n) cost(q_(k)). This also implies that the length of the circuit corresponding to U_(q) can be upper bounded by the function linear in cost vector (L₁, . . . , L_(M)). 3.4. Generalized Lipschitz Order. The Lipschitz order L in the quaternion algebra

$\left( \frac{{- 1},{- 1}}{\mathbb{Q}} \right)$ can be expressed in the following way: L=e ₁ ⁻¹(

[i])+e ₂ ⁻¹(

[i])=

+

i+

j+

k, where

[i] is a maximal order of

(√{square root over (−1)}). For the definition of e₁,e₂ see eq. (4) in the “Using quaternions to represent unitaries” Section (Section 3.1). This construction can be generalized to arbitrary totally definite quaternion algebra. Let

_(K) be a ring of integers of K. It is a two dimensional

_(F) module, therefore it has a

_(F) pseudo basis and can be written (in modified Hermite Normal Form), see Corollary 2.2.9 in Henri Cohen, “Advanced Topics in Computational Number Theory,” Graduate Texts in Mathematics (Springer New York, 2000) as:

_(K)=

_(F) +γI where I is integral

_(F) ideal and γ an element of K such that 1, γ is a F-basis of K.

The generalized Lipschitz order can be defined as: L=e ₁ ⁻¹(

_(K))+e ₂ ⁻¹(

_(K)) L=

_(F) +e ₁ ⁻¹(γ)I+

_(F) j+e ₁ ⁻¹(γ)jI  (12) 3.5. Lattices. In the disclosure, lattices are used that are related a) to the ring of integers

_(K) of CM field K, b) to ideals in

_(K), and c) to the unit group of a the totally real subfield of F of K. Definitions related to integer lattices are briefly reviewed here. More detailed discussion of the definitions and related results can be found in Daniele Micciancio and Sha Goldwasser, “Complexity of lattice problems: a cryptographic perspective,” volume 671 of Springer International Series in Engineering and Computer Science (Springer, 2002).

Let B={b₁, . . . , b_(n)} be a set of linearly-independent vectors in

^(m), where m≥n. The discrete group

(B)=B(

^(n))=

b₁+ . . . +

b_(n) is called the integer lattice of rank n with basis B. Let span (

)=

be the real span of an n-dimensional lattice

and write span(

(B))=span(B). A subset

∈span(

) is called a fundamental domain of the lattice

if for every vector t∈span(

) there exists a unique lattice vector v(t)∈

such that t−v(t)∈

There are at least two different centrally-symmetric fundamental domains associated with each lattice basis. The centered fundamental parallelepiped C(B) associated to a lattice basis B is given by the inequalities C(B)=B[−½,½)^(n) ={Bx:−½≤x _(k)<½,k=1, . . . ,n}. The second fundamental domain is defined in terms of Gram-Schmidt orthogonalization (GSO) B*=[b₁*, . . . , b_(n)*] of a lattice basis B:

$\begin{matrix} {{b_{1}^{*} = b_{1}}{{b_{i}^{*} = {b_{i} - {\sum\limits_{j < i}\;{\mu_{i,j}b_{j}^{*}}}}},}} & (13) \end{matrix}$ where the orthogonalization coefficients μ_(i,j) are defined as μ_(i,j) =

b _(i) ,b _(j) *

/

b _(j) *,b _(j)*

. Note that GSO of a lattice basis is not necessary a basis of

(B). It is related to the original basis via

$B = {{B^{*}\begin{pmatrix} 1 & \; & \mu_{ji} \\ \; & \ddots & \; \\ 0 & \; & 1 \end{pmatrix}}.}$

The centered orthogonalized fundamental parallelepiped C(B*) associated to a lattice basis B is given by the inequalities C(B*)=B*[−½,½)^(n) ={B*x:−½≤x _(k)<½,k=1, . . . ,n} where B* is the Gram-Schmidt orthogonalization of B.

For every integer lattice of rank n>1, there are infinitely many choices of bases. Indeed, for any transformation G∈GL_(n)(

), the basis BG spans the same set of vectors over

as the basis B. Here,

x,y

are used for the standard Euclidean inner product of vectors x,y∈

^(m), and ∥x∥ for the corresponding norm. Reduced lattice bases obtained using the Lenstra-Lenstra-Lovasz (LLL) or Hermite-Korkine-Zolotaroff (HKZ) reduction algorithms allow one to ensure that the sizes of the above fundamental domains are essentially independent of an initial choice of basis, depending on only on the lattice determinant. See, e.g., Arjen K. Lenstra, Hendrik W. Lenstra, and Laszlo Lovasz, “Factoring polynomials with rational coefficients,” Math. Ann., 261:515-534 (1982); Phong Q. Nguyen and Damien Stehle, “Floating-point LLL revisited,” in 24th Annual International Conference on the Theory and Applications of Cryptographic Techniques (EUROCRYPT'05), pages 215-233 (2005).

Elements of

_(K) correspond to 2d dimensional real vectors via map σ:

_(K)→

^(2d) z

(Re(σ_(1,+)(z)),Im(σ_(1,+)(z)), . . . ,Re(σ_(1,+)(z)),Im(σ_(1,+)(z))). The image

of

_(K) under σ is a 2d dimensional integer lattice with associated bilinear form given by Tr_(K/)

(xy*). Each integral basis of

_(K) corresponds to the basis of

. Similarly each

_(K) ideal has

basis. The images of

_(K) ideals under map σ correspond to a sublattices of

. Determinant of

is equal to the discriminant of

_(F).

4. Approximation Algorithm

4.1. High Level Description of the Algorithm.

In this section, a formal description of an exemplary embodiment for solving Problem 2.4 is given and for implementing item 1302 in FIG. 13. A description of the principles underlying the exemplary embodiment are also given. This section is organized as following: each procedure presented in the section is accompanied with the Theorem or Proposition that proves its correctness and bounds on the runtime. In the proofs of the theorems, reference is made to variables defined in the pseudo-code of the corresponding procedures. One of the goals of the remainder of this disclosure is to prove the following theorem:

Theorem 4.1.

There exist constants C_(min),C_(max) and an algorithm (the online part of the procedure APPROXIMATE, FIG. 1) that given

-   -   quaternion gate set specification (see Definition 2.2),     -   real numbers φ and ε∈(0, ½),     -   cost vector (L₁, . . . , L_(M)),         such that     -   narrow class of         ₁ ^(L) ¹ · . . . ·         _(M) ^(L) ^(M) is trivial and     -   L₁ log(N(         ₁))+ . . . +L_(M) log(N(         _(M)))−4 log(1/ε)∈[C_(min),C_(max)]         finds quaternion q from the generalized Lipschitz order such         that nrd(q)         _(F)=         ₁ ^(L) ¹ · . . . ·         _(M) ^(L) ^(M) and d₂(U_(q),R_(z)(φ))≤ε. Constants         C_(min),C_(max) depend only on the quaternion gate set         specification and can be computed in advance (by the offline         part of procedure APPROXIMATE, FIG. 1). The runtime of the         algorithm is on average polynomial in log(1/ε) under the         Conjecture 4.8.

Before looking at details of the algorithm on FIG. 1, the convention used for all pseudo-code shown in the paper is discussed. The example procedures have an offline and an online part. The offline part of the procedures can be executed once (e.g., only once) for a given quaternion gate set specification. Its input is denoted by words Fixed input. For example, for procedure APPROXIMATE the fixed input consists of the most part of quaternion gate set specification (see Definition 2.2) and technical parameter P related to the method for solving relative norm equation. The output of the offline part of the procedure is denoted by term Specification. For example, in the case of procedure APPROXIMATE, the output of the offline part is the additive constant defining the quality of approximation.

The online part of all procedures is executed for each instance of the approximation problem being solved. The instance of the problem is defined by angle φ, target precision ε and target cost vector (L₁, . . . , L_(M)). These are the inputs for the online part of procedure APPROXIMATE. The input to the online part of each procedure is denoted by word Input. An online part of each procedure uses results of computations done in the offline part. Naturally, any offline part can not depend on the results of online computation. The output of online part of each procedure is denoted by word Output. In the complexity analysis discussed here, the primary concern is the online part and it will be shown that the online part of procedure APPROXIMATE has a polynomial runtime on average under a certain number theoretic conjecture. In the “End to end examples of using the framework” Section, the runtime of both online and offline parts is provided for some examples, demonstrating that the offline part is not prohibitively expensive for instances of the problem interesting for applications. Next, the online part of procedure APPROXIMATE is described in more detail.

The method (algorithm) shown on FIG. 1 finds a quaternion q from the generalized Lipschitz order that has two following properties:

-   -   (1) d₂(U_(q),R_(z)(φ))≤ε,     -   (2) nrd(q)         _(F)=         ₁ ^(L) ¹ · . . . ·         _(M) ^(L) ^(M) .         As discussed in the “Generalized Lipschitz order” subsection,         each quaternion q from the generalized Lipschitz order can be         represented using two elements z₁,z₂ of the ring of integers         _(K) of the CM-field K in the following way:         q=e ₁ ⁻¹(z ₁)+e ₂ ⁻¹(z ₂).         As discussed in the “Distance to R_(z) rotations” Section, the         distance d₂(U_(q),R_(z)(φ) depends only on z₁. For this reason,         in the algorithm, z₁ (procedure RANDOM-INTEGER-POINT, line 9 in         FIG. 1) is first picked such that condition (1) above is         satisfied and then z₂ (procedure FAST-SOLVE-NORM-EQ, line 12 in         FIG. 1) is found such that condition (2) is also satisfied.

Procedure FAST-SOLVE-NORM-EQ solves the relative norm equation N_(K/F)(z₂)=z₂z₂=e in the relative extension K/F for a special class of right hand sides e. There are two challenges related to this procedure. First, the solution does not always exist for arbitrary right hand side e. Second, solving the arbitrary instance of the norm equation (for fixed extension K/F) can be as hard as factoring. Both these challenges can be addressed. First, the necessary conditions on the right hand side e of the equation to be solvable are identified. Second, the set of right hand sides e for which the equation can be solved in probabilistic polynomial time (using procedure IS-EASILY-SOLVABLE, line 11, in FIG. 1) is identified and attempt(s) are made to solve the equation only for such right hand sides. The claim that example embodiments of the disclosed technology work in polynomial time is conditional on the conjecture that the procedure IS-EASILY-SOLVABLE returns true with probability O(1/log(1/ε)). This conjecture is discussed in more detail in the “On the conjecture related to the approximation algorithm performance” Section (Section 4.5). For more details, see also the “Relative norm equations” Section (Section 5).

The necessary condition for the norm equation to be solvable is that for all embeddings σ_(k) of F into

it must be the case that σ_(k)(e)>0 (for k=1, . . . , d). Procedure RANDOM-INTEGER-POINT (line 9 in FIG. 1) returns z₁ such that these conditions are satisfied together with condition d₂(U_(q),R_(z)(φ))≤ε. Recall, that one can associate an integer lattice

with

_(K) (see the “Lattices” Section (Section 3.5)). Essentially, procedure RANDOM-INTEGER-POINT samples lattice points from the following convex subset of

^(2d): S _(r,φε) ={x∈

^(2d):Re((x ₁ +ix ₂ −z ₀ ^(r,φ,ε))e ^(−iφ/2))≥0,|x _(2k) +ix _(2k+1)|²≤σ_(k)(r)}  (14) where z₀ ^(r,φ,ε)=√{square root over (σ₁(r))}(1−ε²)e^(−iφ/2). See FIG. 14 for the visualization of the projection of the set above. The geometry of the set S_(r,φ,ε) is determined by the reduced norm r of the quaternion output by procedure APPROXIMATE. One has freedom in choosing r up to a unit of

_(F). It can then be used to find r such that one can sample lattice points from S_(r,φ,ε) in time that is logarithmic in the volume of S_(r,φ,ε). Such r is found by procedure SUITABLE-Q-NORM (line 6 in FIG. 1).

Before proceeding to the proof of Theorem 4.1, results that are proven later and used in the proof are presented.

Theorem 4.2.

There exist real numbers p₀, M and vectors R^(min),R^(max) from (0, ∞)^(d) (computed in the offline part of RANDOM-INTEGER-POINT procedure, FIG. 2) such that for any real number φ, real number ε∈(0, ½) and totally positive r from

_(F) such that: √{square root over (σ₁(r))}ε²/4∈[R ₁ ^(min) ,R ₁ ^(max)], √{square root over (σ_(k)(r))}∈[R _(k) ^(min) ,R _(k) ^(max)] for k=2, . . . ,d there is an element z from

_(F) such that σ(z) is in S_(r,φ,ε) (see Equation (14)). Procedure RANDOM-INTEGER-POINT runs in time polynomial in log(1/ε) and returns each element of the set Cand _(r,φ,ε) ={z∈

_(K):σ(z)∈S _(r,φ,ε)}  (15) with probability at least p₀/|Cand_(r,φ,ε)|. The size of the set Cand_(r,φ,ε) belongs to the interval [2√{square root over (4−ε²)}/ε,(4√{square root over (4−4ε²)}/ε+2)M]

The proof of Theorem 4.2 can be found in the “Solution region sampling” Section (Section 4.3).

Theorem 4.3.

Given vector R^(min) from (0, ∞)^(d), there exists constants C_(min),C_(max) and vector R^(max) (computed by the offline part of SUITABLE-Q-NORM procedure, FIG. 3) such that for all non-negative integers L₁, . . . , L_(M) and a positive real number ε∈(0, 1) satisfying: L ₁ log N(

₁)+ . . . +L _(M) log N(

_(M))−4 log(1/ε)∈[C _(min) ,C _(max)]

there is an algorithm that decides if narrow class number of ideal

₁ ^(L) ¹ · . . . ·

_(M) ^(L) ^(M) is trivial. If this is the case, the algorithm outputs totally positive element r of

_(F) such that r

_(F)=

₁ ^(L) ¹ · . . . ·

_(M) ^(L) ^(M) and √{square root over (σ₁(r))}ε²/4∈[R ₁ ^(min) ,R ₁ ^(max)], √{square root over (σ_(k)(r))}∈[R _(k) ^(min) ,R _(k) ^(max)] for k=2, . . . ,d The algorithm runs in polynomial time in log(1/ε) (see online part of SUITABLE-Q-NORM procedure in FIG. 3).

The informal discussion and the proof of the Theorem 4.3 can be found in the “Picking a suitable quaternion norm” Section (Section 4.2).

Theorem 4.4.

Given totally positive element e of F, there exists an algorithm for testing if the instance of integral relative norm equation in K/F zz*=e,z∈

_(K) can be solved in polynomial time in log T₂(e) (procedure IS-EASILY-SOLVABLE, where T₂(e)=Σ_(k=1) ^(d)σ_(k) ²(e)). If the test is passed, there exist another algorithm for deciding if the solution exists and finding it that runs in time polynomial in log T₂(e) (procedure FAST-SOLVE-NORM-EQ). Procedure IS-EASILY-SOLVABLE returns true for at least those cases when the ideal e

_(F) is prime. Proof of Theorem 4.1. First, a proof is given that the output is correct. The norm of the quaternion q computed on line 15 (FIG. 1) is equal to totally positive element r of

_(F) computed on line 6 (FIG. 1). This is because on line 12 (FIG. 1) one finds z₂ such that r=nrd(q)=|z₁|²−b|z₂|². By Theorem 4.3 the output r of the procedure SUITABLE-Q-NORM (line 6, FIG. 1) satisfies r

_(F)=nrd(q)

_(F)=

₁ ^(L) ¹ · . . . ·

_(M) ^(L) ^(M) . According to Theorem 4.2 procedure RANDOM-INTEGER-POINT (line 9, FIG. 1) returns an element z₁ of

_(K) such that Re((σ_(1,+)(z)−z ₀)e ^(−iφ/2))≥0,|σ_(1,+)(z)|≤√{square root over (σ₁(r))} where z₀=√{square root over (σ₁(r))}(1−ε²)e^(−iφ/2). According to Proposition 3.1 this implies that d₂(R_(z)(φ),U_(q))≤ε. Numbers z₁,z₂ are in

_(F) which immediately implies that q is in generalized Lipschitz order.

Next, it is shown that the restrictions on inputs of all the procedures called within online part of procedure APPROXIMATE are satisfied. Procedure SUITABLE-Q-NORM (line 6, FIG. 1) always succeeds if the narrow class group of

₁ ^(L) ¹ · . . . ·

_(M) ^(L) ^(M) is trivial and L ₁ log N(

₁)+ . . . +L _(M) log N(

_(M))−4 log(1/ε)∈[C _(min) ,C _(max)], which is required in the statement of the theorem. From Theorem 4.3 it is known that procedure SUITABLE-Q-NORM (line 6, FIG. 1) finds r such that: √{square root over (σ₁(r))}ε²/4∈[R ₁ ^(min) ,R ₁ ^(max)], √{square root over (σ_(k)(r))}∈[R _(k) ^(min) ,R _(k) ^(max)] for k=2, . . . ,d therefore procedure RANDOM-INTEGER-POINT always succeeds. Now it is shown that e=(r−z₁z₁*)/(−b) is totally positive, where b is the parameter from the definition of the quaternion algebra Q. Note that e being totally positive is required by procedures IS-EASILY-SOLVABLE (line 11, FIG. 1) and FAST-SOLVE-NORM-EQUATION (line 12, FIG. 1). For all k=1, . . . , d it the case that σ_(k)(e)=(σ_(k)(r)−|σ_(k,+)(z ₁)|²)/−σ_(k)(b). By definition, b is totally negative and σ_(k)(e)>0 if and only if σ_(k)(r)−|σ_(k,+)(z₁)|²>0. By Theorem 4.2, the output of procedure RANDOM-INTEGER-POINT (line 9, FIG. 1) satisfies |σ_(k,+)(z₁)|≤R_(k) ^(min) for k=2, . . . , d. By Theorem 4.3 totally positive element r of

_(F) satisfies R_(k) ^(min)≤√{square root over (σ_(k)(r))} for k=2, . . . , d and therefore σ_(k)(e)>0 for k=2, . . . , d. It was already shown above that |σ_(k,+)(z₁)|≤√{square root over (σ₁(r))} which implies σ₁(e)>0.

It remains to show that the example algorithm terminates and runs on average in time polynomial in log(1/ε). Procedure SUITABLE-Q-NORM (line 6, FIG. 1) runs in time polynomial in log(1/ε) by Theorem 4.3. Now it is shown that all procedures inside the loop run in polynomial time. Procedure RANDOM-INTEGER-POINT runs in polynomial time in log(1/ε) according to Theorem 4.2.

Next, it is shown that the logarithm of

${T_{2}(e)} = {\sum\limits_{k = 1}^{d}{\sigma_{k}^{2}(e)}}$ is bounded by polynomial in log(1/ε). This implies that procedures IS-EASILY-SOLVABLE and FAST-SOLVE-NORM-EQ run on average in polynomial time according to Theorem 4.4. Indeed, one has σ_(k)(e)≤σ_(k)(r)/σ_(k)(−b), √{square root over (σ₁(r))} is bounded by 4R₁ ^(max)/ε² and √{square root over (σ_(k)(r))} are bounded by R_(k) ^(max) for k=2, . . . , d. Finally, arithmetic in the number field (line 10, FIG. 1) can be done in time polynomial in log T₂(z₁) and log T₂ (r)·log T₂(z₁) is also bound using inequalities |σ_(k+1)(z₁)|²≤σ_(k)(r). Also performed is the test e∈

_(F) (line 11, FIG. 1). It can be performed in polynomial time in log T₂(r), because it can be reduced to multiplying a vector over

with the norm bounded by C log T₂(r) by fixed matrix over

.

It can be concluded that procedure APPROXIMATE runs on average in time polynomial in log(1/ε) under the conjecture that the fraction of points in set {z∈

_(K):σ(z)∈S_(r,φ,ε)} for which one can reach line 12 (FIG. 1) and successfully find z₂ scales as Ω(log(1/ε)). This is because according to Theorem 4.2 all points from {z∈

_(K):σ(z)∈S_(r,φ,ε)} are sampled sufficiently uniformly.

4.2. Picking a Suitable Quaternion Norm.

In this subsection, the following theorem is proven:

Theorem 4.3.

Given vector R^(min) from (0, ∞)^(d), there exists constants C_(min),C_(max) and vector R^(max) (computed by the offline part of SUITABLE-Q-NORM procedure, FIG. 3) such that for all non-negative integers L₁, . . . , L_(M) and a positive real number ε∈(0, 1) satisfying: L ₁ log N(

₁)+ . . . +L _(M) log N(

_(M))−4 log(1/ε)∈[C _(min) ,C _(max)] there is an algorithm that decides if narrow class number of ideal

₁ ^(L) ¹ · . . . ·

_(M) ^(L) ^(M) is trivial. If this is the case, the algorithm outputs totally positive element r of

_(F) such that r

_(F)=

₁ ^(L) ¹ · . . . ·

_(M) ^(L) ^(M) and √{square root over (σ₁(r))}ε²/4∈[R ₁ ^(min) ,R ₁ ^(max)], √{square root over (σ_(k)(r))}∈[R _(k) ^(min) ,R _(k) ^(max)] for k=2, . . . ,d The algorithm runs in polynomial time in log(1/ε) (see online part of SUITABLE-Q-NORM procedure in FIG. 3).

The proof relies on the following proposition proven in this and the next sections.

Proposition 4.5.

Given non-negative integers L₁, . . . , L_(M) there is an algorithm (procedure TOTALLY-POS-GEN, FIG. 7) that decides if ideal

₁ ^(L) ¹ · . . . ·

_(M) ^(L) ^(M) has a totally positive generator. The algorithm also finds a totally positive generator r of the ideal if it exists. The algorithm runs in time polynomial in L₁, . . . , L_(M) and log T₂(r) is bounded by the function that is linear in L₁, . . . , L_(M). Proposition 4.6. There exists real numbers δ₀ and δ₁, . . . , δ_(d) (computed by the offline part of the procedure UNIT-ADJUST, FIG. 9) such that there exists an algorithm (online part of the procedure UNIT-ADJUST, FIG. 9) that for any real numbers t₁, . . . , t_(d) finds a unit u of

_(F) such that the following inequalities hold |log|σ_(k)(u)|−t _(k)|≤log δ_(k) for k=1, . . . ,d, under the assumption that |t₁+ . . . +t_(d)<log δ₀. The runtime of the algorithm is bounded by a polynomial in ∥t∥. Proposition 4.7. Given real numbers δ₀, . . . , δ_(d)>1, vector R^(min) from (0, ∞)^(d), and prime ideals

₁, . . . ,

_(M) there exist real numbers C_(min),C_(max) and vector R^(max) from (0, ∞)^(d) (computed by the offline part of the procedure TARGET-SIZE, FIG. 8) such that there exist an algorithm (online part of the procedure TARGET-SIZE, FIG. 8) that given non-negative integers L₁, . . . , L_(M), real number ε and totally positive element r of F finds real numbers t₁, . . . , t_(d) such that t ₁+log(√{square root over (σ_(k)(r))}ε²/4)∈[log R ₁ ^(min)+log δ₁,log R ₁ ^(max)−log δ₁] t _(k)+log √{square root over (σ_(k)(r))}∈[log R _(k) ^(min)+log δ_(k),log R _(k) ^(max)−log δ_(k)], k=2, . . . ,d and |t₁+ . . . +t_(d)|<log δ₀. The algorithm succeeds under the assumption that L₁ log N(

₁)+ . . . +L_(M) log N(

_(M))−4 log(1/ε)∈[C_(min),C_(max)].

The runtime of the algorithm is bounded by a polynomial in log(1/ε) and log T₂(r). The norm ∥t∥ is bounded by the function that is linear in the same variables.

Proof of Theorem 4.3.

First, it is proven that the procedure terminates in polynomial time in L₁, . . . , L_(M) when it returns FALSE. Indeed, when ideal

₁ ^(L) ¹ · . . . ·

_(M) ^(L) ^(M) does not have a totally positive generator procedure TOTALLY-POS-GENERATOR returns FALSE and procedure SUITABLE-Q-NORM terminates (line 8, FIG. 3). According to Proposition 4.5 procedure TOTALLY-POS-GENERATOR runs in polynomial time in L₁, . . . , L_(M).

Next, the case when the output of SUITABLE-Q-NORM procedure is TRUE is considered. First, it is proven that the output of SUITABLE-Q-NORM procedure (FIG. 3) is correct in this case. By Proposition 4.5, algebraic integer r (line 6, FIG. 3) is totally positive and r

_(F)=

₁ ^(L) ¹ · . . . ·

_(M) ^(L) ^(M) . By Proposition 4.6, algebraic integer u (line 11, FIG. 3) is a unit. Therefore ru² is also totally positive and also generates ideal

₁ ^(L) ¹ · . . . ·

_(M) ^(L) ^(M) .

It remains to show that σ_(k)(ru²) satisfy required inequalities. By Proposition 4.6 unit u (line 11, FIG. 3) computed by procedure UNIT-ADJUST satisfies the following inequalities: |log|σ_(k)(u)|−t _(k)|≤log δ_(k) , k=1, . . . ,d, because by Proposition 4.7 procedure TARGET-SIZE ensures that |t₁+ . . . +t_(d)|<log δ₀. Now, it can be seen that log √{square root over (σ_(k)(ru ²))}−log √{square root over (σ_(k)(r))}=log|σ_(k)(u)|∈[t _(k)−log δ_(k) ,t _(k)+log δ_(k)] log(√{square root over (σ_(k)(ru ²))}ε²/4)−log(√{square root over (σ₁(r))}ε²/4)=log|σ₁(u)|∈[t ₁−log δ₁ ,t ₁+log δ₁] This immediately implies that log √{square root over (σ_(k)(ru ²))}∈[log R _(k) ^(min),log R _(k) ^(max)], for k=2, . . . ,d log(√{square root over (σ₁(ru ²))}ε²/4)∈[log R ₁ ^(min),log R ₁ ^(max)] It is now shown that the runtime of the algorithm is bounded by a polynomial in log(1/ε). All L_(k) are bounded by function linear in log(1/ε). Procedure TOTALLY-POS-GENERATOR runs in polynomial time and produces r such that log T₂ (r) is bounded by a function linear in L₁, . . . , L_(M). This ensures that procedure TARGET-SIZE outputs t₁, . . . , t_(d) such that their bit size is bounded by polynomial in log(1/ε). It also ensures that ∥t∥ is bounded by a function linear in log(1/ε). This ensures that procedure UNIT-ADJUST runs in polynomial time. Note that for unit u (computed in line 11, FIG. 3) one has

${\log\sqrt{T_{2}(u)}} = {{\log\sqrt{\sum\limits_{k = 1}^{d}{\sigma_{k}^{2}(u)}}} \leq {\sqrt{d}{\max\limits_{k}{\log{\sigma_{k}}}}} \leq {\sqrt{d}{\max\limits_{k}\left( {t_{k} + {\log\mspace{11mu}\delta_{k}}} \right)}} \leq {\sqrt{d}\left( {{t} + {\max\limits_{k}{\log\mspace{11mu}{\delta_{k}.}}}} \right)}}$ Therefore log T₂(u) is bounded by a function linear in log(1/ε). Hence, the time spent on computing ru² is bounded by polynomial in log(1/ε). It has therefore been shown that procedure SUITABLE-Q-NORM runs in polynomial time.

Now, the procedure TOTALLY-POS-GENERATOR (FIG. 7) is discussed in more detail. Finding totally positive generator of the ideal is strictly more difficult than finding a generator of the ideal. The latter problem is known to be hard and there is no polynomial time algorithm known for it. Here, because the ideal of interest has special form

₁ ^(L) ¹ . . .

_(m) ^(L) ^(M) one can find a totally positive generator for it by precomputing certain information about

₁, . . . ,

_(m) (in the offline part of TOTALLY-POS-GENERATOR procedure, FIG. 7). Now, the following proposition is proven:

Proposition 4.5.

Given non-negative integers L₁, . . . , L_(M) there is an algorithm (procedure TOTALLY-POS-GEN, FIG. 7) that decides if ideal

₁ ^(L) ¹ · . . . ·

_(M) ^(L) ^(M) has a totally positive generator. The algorithm also finds a totally positive generator r of the ideal if it exists. The algorithm runs in time polynomial in L₁, . . . , L_(M) and log T₂(r) is bounded by the function that is linear in L₁, . . . , L_(M). Proof. First, the correctness of the online part of the algorithm is proven. It is not difficult to see that

₁ ^(L) ¹ · . . . ·

_(M) ^(L) ^(M) =(

₁ ^(N) ¹ )^(s) ¹ · . . . ·(

_(M) ^(N) ^(M) )^(s) ^(M)

₁ ^(l) ¹ · . . . ·

_(M) ^(l) ^(M) One knows what each of ideals

_(k) ^(N) ^(k) has totally positive generator. Therefore

₁ ^(L) ¹ · . . . ·

_(M) ^(L) ^(M) has totally positive generator if and only if ideal

₁ ^(l) ¹ · . . . ·

_(M) ^(l) ^(M) does. This is what is checked on line 11 in FIG. 7. It has been shown that the procedure always returns FALSE, when

₁ ^(L) ¹ · . . . ·

_(M) ^(L) ^(M) does not have totally positive generator, and TRUE when it does. In the case when the procedure returns TRUE, r[l₁, . . . , l_(k)]

_(F)=

₁ ^(l) ¹ · . . . ·

_(M) ^(l) ^(M) and therefore

₁ ^(L) ¹ · . . . ·

_(M) ^(L) ^(M) =r ₁ ^(s) ¹

_(F) · . . . ·r _(M) ^(s) ^(M)

_(F) ·r[l ₁ , . . . ,l _(k)]

_(F) It can be seen that r₁ ^(s) ¹ . . . r_(M) ^(s) ^(M) r[l₁, . . . , l_(k)] is a totally positive element of

_(F) that generates ideal

₁ ^(L) ¹ · . . . ·

_(M) ^(L) ^(M) .

Now, it is shown that the algorithm runs in polynomial time. The number of multiplications to perform is bounded by L₁+ . . . +L_(M). Note that

${{T_{2}\left( {x_{1} \cdot \ldots \cdot x_{n}} \right)} \leq {d{\prod\limits_{k = 1}^{n}{T_{2}\left( x_{k} \right)}}}},{{{for}\mspace{14mu} x_{k}} \in {F.}}$ Therefore, each time a multiplication is performed, the value log T₂ of the arguments is bounded by a function linear in L₁, . . . , L_(M). It can be concluded that r₁ ^(s) ¹ . . . r_(M) ^(s) ^(M) r[l₁, . . . , l_(k)] is computed in time polynomial in L₁, . . . , L_(M). The inequality above also implies that log T₂ of the algorithm is bounded by a function linear in L₁, . . . , L_(M).

It remains to show the correctness of the offline part of procedure. First, note that N_(k) always exist. The fact that the class group of the number field is always finite implies that for each ideal

_(k) there exists a number N_(k)′ (dividing the order of class group) such that ideal

_(k) ^(N′) ^(k) is principal and has generator ξ_(k). Note that ξ_(k) ² is totally positive generator of

_(k) ^(2N′) ^(k) . It has been shown that N_(k)≤2N_(k)′. Computing the order of the class group and testing if the ideal is principal are two standard problems solved in computational number theory. Checking that given ideal has a totally positive generator can also be done using standard methods. If totally positive generator r exists then the ratio of the ideal generator ξ and r must be a unit. Therefore, once ξ is known it remains to find unit u such that ξu is totally positive. This can be done if and only if ξ

_(F) has a totally positive generator. It is not difficult to check this algorithmically, once the system of fundamental units of

_(F) has been computed.

Note, that in case of Clifford+T and Clifford+R_(z)(π/16) gate set the narrow class group of F is trivial and therefore any ideal has a totally positive generator. This significantly simplifies procedure TOTALLY-POS-GENERATOR.

On a high level, procedure TARGET-SIZE (FIG. 8) is about finding a solution to the system of linear inequalities for t₁, . . . , t_(d) mentioned in the statement of the Proposition 4.7. In the offline part, constant C that guarantees that the system has a solution is computed. In addition, it will be shown that the solution to such a system can be represented by numbers of moderate bit-size, can be found in polynomial time and that the norm of the solution vector is bounded. Now, the formal proof is presented.

Proposition 4.7.

Given real numbers δ₀, . . . , δ_(d)>1, vector R^(min) from (0, ∞)^(d), and prime ideals

₁, . . . ,

_(M) there exist real numbers C_(min),C_(max) and vector R^(max) from (0, ∞)^(d) (computed by the offline part of the procedure TARGET-SIZE, FIG. 8) such that there exist an algorithm (online part of the procedure TARGET-SIZE, FIG. 8) that given non-negative integers L₁, . . . , L_(M), real number ε and totally positive element r of F finds real numbers t₁, . . . , t_(d) such that t ₁+log(√{square root over (σ_(k)(r))}ε²/4∈[log R ₁ ^(min)+log δ₁,log R ₁ ^(max)−log δ₁] t _(k)+log √{square root over (σk(r))}∈[log R _(k) ^(min)+log δ_(k),log R _(k) ^(max)−log δ_(k)], k=2, . . . ,d and |t₁+ . . . +t_(d)|<log δ₀. The algorithm succeeds under the assumption that L₁ log N(

₁)+ . . . +L_(M) log N(

_(M))−4 log(1/ε)∈[C_(min),C_(max)].

The runtime of the algorithm is bounded by a polynomial in log(1/ε) and log T₂(r). The norm ∥t∥ is bounded by the function that is linear in the same variables.

Proof.

First, the correctness of the procedure TARGET-SIZE (FIG. 8) is proven. On lines 6, 7 in FIG. 8, the smallest values for which the constraints on t_(k) are satisfied are assigned to t_(k) Next, it is shown that the choice of C_(min), C_(max) and R^(max) ensures that all other constraints are also satisfied. Indeed, R_(k) ^(max)≥R_(k) ^(min)+2 log δ_(k) therefore intervals for t_(k) are non-empty. It remains to show that |t₁+ . . . +t_(d)|≤δ₀. This follows from the following equalities:

$\begin{matrix} {{t_{1} + \ldots + t_{d}} = {{\sum\limits_{k = 1}^{d}{\log\mspace{11mu} R_{k}^{\min}}} + {\sum\limits_{k = 1}^{d}{\log\mspace{11mu}\delta_{k}}} - {\log\left( {\sqrt{N_{F/{\mathbb{Q}}}(r)}{ɛ^{2}/4}} \right)}}} \\ {= {{C_{\min}/2} - {\frac{1}{2}{\sum\limits_{k = 1}^{M}{L_{k}\mspace{11mu}\log\mspace{11mu}{N\left( {\mathfrak{p}}_{k} \right)}}}} + {2\mspace{11mu}{\log\left( {1/ɛ} \right)}}}} \end{matrix}$ Therefore value t₁+ . . . +t_(d) belongs to the interval [(C_(min)−C_(max))/2, 0]. The definition of C_(max) precisely implies that [(C_(min)−C_(max))/2, 0]=[−log δ₀, 0] which gives required bound on the sum of t_(k). Note that the analysis above performed for t_(k) will hold for any t_(k)′ in the interval [t_(k),t_(k)+log δ₀/d].

Parameter log δ₀ is needed to account for the finite precision arithmetic used. It is not difficult to see that as soon as precision of arithmetic used is smaller then log δ₀/C₁d for sufficiently big fixed constant C₁ numbers t_(k) computed within mentioned precision will satisfy all required constraints. It is sufficient to perform the calculation up to fixed precision independent on the online part of the algorithm input. This implies that all calculations in the online part can be performed in polynomial time. Also note that t_(k) are bounded by functions linear in log T₂(r) and log(1/ε) and therefore the same is true for ∥t∥. As it has been established that t_(k) can be computed up to fixed precision, bound on ∥t∥ implies a bound on the number of bits needed to specify each t_(k). This concludes the proof.

4.3. Solution Region Sampling.

In this subsection, the Theorems 4.2 used in the proof of Theorem 4.1 in Section 4.1 are proven:

Theorem 4.2.

There exist real numbers p₀, M and vectors R^(min),R^(max) from (0, ∞)^(d) (computed in the offline part of RANDOM-INTEGER-POINT procedure, FIG. 2) such that for any real number φ, real number ε∈(0, ½) and totally positive r from

_(F) such that: √{square root over (σ₁(r))}ε²/4∈[R ₁ ^(min) ,R ₁ ^(max)], √{square root over (σ_(k)(r))}∈[R _(k) ^(min) ,R _(k) ^(max)] for k=2, . . . ,d there is an element z from

_(F) such that σ(z) is in S_(r,φ,ε) (see Equation (14)). Procedure RANDOM-INTEGER-POINT runs in time polynomial in log (1/ε) and returns each element of the set Cand _(r,φ,ε) ={z∈

_(K):σ(z)∈S _(r,φ,ε)}  (15) with probability at least p₀/|Cand_(r,φ,ε)|. The size of the set Cand_(r,φ,ε) belongs to the interval [2√{square root over (4−ε²)}/ε,(4√{square root over (4−4ε²)}/ε+2)M] Proof. The procedure returns points from the set Cand_(r,φ,ε) because on line 21 in FIG. 2, z is selected such that σ(z)∈S_(r,φ,ε) (see Equation (15)). It is now shown that the procedure terminates, on average, after a fixed number of loop iterations. The probability p_(Cand) of variable Sample-found being TRUE ca be estimated (see line 21 in FIG. 2). The number of the main loop iterations of RANDOM-INTEGER-POINT procedure has a geometric distribution with parameter p_(Cand). Now consider lower bound p_(Cand). Let H′=Rε√{square root over (4−ε²)} (see FIG. C). When the absolute value of N (on line 16 in FIG. 2) is less then ┌H′/(ε²R)┐−1 (see FIG. 15) the shifted fundamental region t+C(B) is a subset of S_(r,φ,ε) (note that t is computed on line 17 in FIG. 2 based on N). In the case if z′=0 (line 19 in FIG. 2) inclusion t+C(B)⊂S_(r,φ,ε) implies that σ(z) is in S_(r,φ,ε). The probability of z′=0 is equal to 1/|SHIFTS| (see line 6 in FIG. 2). For this reason one can lower bound p_(Cand) as:

$p_{Cand} \geq {\frac{{2\left\lceil {H^{\prime}/\left( {ɛ^{2}R} \right)} \right\rfloor} - 1}{{2\left\lceil {H/\left( {ɛ^{2}R} \right)} \right\rfloor} + 1} \cdot \frac{1}{{SHIFTS}}} \geq {\frac{{H^{\prime}/H} - {2\left( {ɛ^{2}{R/2}H} \right)}}{1 + {2\left( {ɛ^{2}{R/2}H} \right)}} \cdot \frac{1}{{SHIFTS}}}$ It is noted that H′/H=√{square root over (4−ε²)}/(2√{square root over (4−4ε²)})∈(½, 1/√{square root over (3)}) and (ε²R/2H′)=/(4√{square root over (4−4ε²)}) (0, 1/(8√{square root over (3)})). The constraint ε∈(0, ½) implies that P_(Cand) is lower bounded by constant independent on the input to the online part of RANDOM-INTEGER-POINT procedure.

FIG. 15 is an illustration 1500 of the sampling scheme. The vertical axis is Im(z) and the horizontal axis is Re(z). Small blue circles (e.g., circles 1510, 1512, 1514) have radius R₁ ^(min). During the example sampling procedure, the center of a possible blue circle is picked at random. Then the sampling procedure finds point Σ_(k)m_(k)z_(k) (see FIG. 2) from

_(K) such that σ_(1,+)(Σ_(k)m_(k)z_(k)) is in blue circle. Next sampling procedure computes z by adding randomly picked shift z′ to Σ_(k)m_(k)z_(k) and checks if σ(z) is in S_(r,φ,ε). The set P₁S_(r,φ,ε) is a blue circle segment 1520 as illustrated in FIG. 15.

It has been shown that the online part of the example procedure comprises a fixed number of arithmetic operations on average. To show that the procedure runs in polynomial time, it is sufficient to show that the absolute value of the logarithm of absolute precision required for the computation is bounded by a polynomial in log(1/ε). Consider line 18 in FIG. 2 in more detail. Let (B⁻¹)′ be an approximation to B⁻¹ to within precision δ_(c) and let t′ be an approximation to t within precision δ_(c) and let m′=(B⁻¹)′t′. Now it will be shown that the norm of the projection of Bm′−t′ can also be bounded in terms of R_(k) ^(min). For convenience, let P_(k) be a projector on the subspace spanned by e_(2k−1),e_(2k). ∥P _(k)(Bm′−t)∥≤∥t−t′∥+∥P _(k) B(┌B ⁻¹)′t′┘−(B ⁻¹)′t′)∥+∥t′∥∥B(B ⁻¹)′−I∥ Now it can be seen that ∥t−t′∥ is bounded by δ_(c), the second term in the sum above is bounded by R_(k) ^(min) and the third term is bounded by some fixed constant times δ_(c)∥t′∥. This implies that one can find m′ such that ∥P_(k)(Bm′−t)∥≤R_(k) ^(min)+δ_(c)′. The absolute value of the logarithm of absolute precision required for the computation is bounded by polynomial in log(1/ε) and log(1/δ_(c)′) because log∥t′∥ is bounded by polynomial in log(1/ε). This is sufficient for purposes of this discussion because it is sufficient to choose log(1/δ_(c)′) to be of order log∥t′∥.

It remains to show that the procedure can get every point from Cand_(r,φ,ε) with probability at least p₀/|Cand_(r,φ,ε)|. First, some notation is introduced for the proof. Let ξ_(N) be a random variable corresponding to N (line 16 in FIG. 2), ξ_(z′) be a random variable corresponding z′ (line 19 in FIG. 2) and ξ_(z) be a random variable corresponding to value z computed on line 20 in FIG. 2. Random variable ξ_(z) is a function of ξ_(N) and ξ_(z′). Random variable ξ_(N) takes integer values in the range [−N_(max),N_(max)] with equal probability and ξ_(z′) takes values in the finite set SHIFTS also with equal probability. Variables ξ_(z′) and ξ_(N) are independent. Let z₀ be a fixed element of Cand_(r,φ,ε). It is sufficient to lower bound P(ξ_(z)=z₀) because the probability of the output of RANDOM-REGION-CENTER being z₀ is P(ξ_(z)=z₀)/P(ξ∈Cand_(r,φ,ε)). Given z₀ there are unique values of N (line 16), m (line 18) and t (line 17), denoted as N₀, m⁰ and t₀, such that σ(z₀) is in S _(r,φ,ε) ∩{x∈

^(2d) :

x−t ₀ ,ΔZ

∈(−½,½]}. The equality

${P\left( {\xi_{z} = z_{0}} \right)} = {{{P\left( {\xi_{z} = {{z_{0}/\xi_{N}} = N_{0}}} \right)}{P\left( {\xi_{N} = N_{0}} \right)}} = \frac{P\left( {\xi_{z} = {{z_{0}/\xi_{N}} = N_{0}}} \right)}{{2N_{\max}} + 1}}$ implies that it is sufficient to lower bound P(ξz=z₀/ξ_(N)=N₀) and relate 2N_(max)+1 to the size of Cand_(r,φ,ε). Note that P(ξ_(z)=z₀/ξ_(N)=N₀)=P(z₀=ξ_(z)′+Σ_(k) m_(k) ⁰z_(k)). This implies that P(ξ_(z)=z₀/ξ_(N)=N₀) is 1/|SHIFTS| if z₀−Σ_(k)m_(k) ⁰z_(k) belongs to set SHIFTS. It is now shown that z₀−ρ_(k)m_(k)z_(k) is always in SHIFTS. It is sufficient to show that

${{P_{k}\left( {{\sigma\left( z_{0} \right)} - {\sigma\left( {\sum\limits_{k}{m_{k}^{0}z_{k}}} \right)}} \right)}} = {{{P_{k}\left( {{\sigma\left( z_{0} \right)} - {Bm}^{0}} \right)}} \leq {R_{k}^{\max} + {R_{k}^{\min}\mspace{14mu}{for}}}}$ k = 2, …  , d ${{P_{1}\left( {{\sigma\left( z_{0} \right)} - {Bm}^{0}} \right)}} \leq {2\sqrt{5}R_{1}^{\max}}$ It is useful to note that ∥P _(k)(σ(z ₀)−Bm ⁰)∥≤∥P _(k)(σ(z ₀)−t ₀)∥+∥P _(k)(t ₀ −Bm ⁰)≤∥P _(k)(σ(z ₀)−t ₀)∥+R _(k) ^(min) The fact that σ(z₀) is in S_(r,φ,ε)∩{x∈

^(2d):

ΔZ

∈(−½,½]} implies that ∥P _(k)(σ(z ₀)−t ₀)∥≤√{square root over (σ_(k)(r))}≤R _(k) ^(max) for k=2, . . . ,d. To establish bound on ∥P₁(σ(z₀)−Bm⁰)∥, observe that P₁σ(z₀) and P₁Bm⁰ both belong to a set with the diameter √{square root over (5)}ε²R/2. It has also been shown that: Cand _(r,φ,ε)⊂SHIFTS+{σ⁻¹(B┌B ⁻¹(Z _(c) +NΔZ)┘):N∈[−N _(max) ,N _(max)]}

Finally, note that: {z∈Z _(K):σ(z)∈Z _(c) +NΔZ+C(B),|N|≤┌H′/(ε² R)┘−1,N∈

}⊂Cand _(r,φ,ε). This implies that if Cand_(r,φ,ε) is non-empty, one has:

$\frac{{Cand}_{r,\varphi,ɛ}}{{2N_{\max}} + 1} \geq \frac{{Cand}_{r,\varphi,ɛ}}{{2{H/\left( {ɛ^{2}R} \right)}} + 2} \geq \frac{{Cand}_{r,\varphi,ɛ}}{{\frac{H}{H^{\prime}}\left( {{2\left\lceil {H^{\prime}/\left( {ɛ^{2}R} \right)} \right\rfloor} + 1} \right)} + 2} \geq \frac{{Cand}_{r,\varphi,ɛ}}{6 + {2{{Cand}_{r,\varphi,ɛ}}}} \geq \frac{1}{8}$ It can be concluded that

${{P\left( {\xi_{z} = z_{0}} \right)} \geq {\frac{1}{{Cand}_{r,\varphi,ɛ}} \cdot p_{0}}},$ where P₀ is

$\frac{1}{8{{SHIFTS}}}.$ Above derivation also gives the required bounds on the size of Cand_(r,φ,ε). 4.3.1. Implementation Aspects. In practice, example implementations of the disclosed technology are looking for the best possible value (or other favorable values) of the additive constants C_(min),C_(max) in Theorem 4.1 while maintaining the polynomial runtime of the online part of the algorithm. Section 8 explains in detail the version of the procedure used in the example implementations used herein to obtain the numerical results reported below (Section 6). In practice, the Nearest Plane Algorithm can be used. See Laszlo Babai, “On Lovasz lattice reduction and the nearest lattice point problem,” Combinatiorica, 6(1):1-13 (1986). It is also possible to show that R_(k) ^(min) can be chosen to be based on C(B*), not based on C(B). One can ensure that the basis used is Hermite-Korkine-Zolotarev reduced which makes it possible to guarantee that R_(k) ^(min) are bounded by some functions of discriminant of

_(K) and these bounds are independent of the choice of the basis of

_(K). See Guillaume Hanrot and Damien Stehle, “Improved analysis of kannan's shortest lattice vector algorithm” in Alfred Menezes, editor, Advances in Cryptology—CRYPTO 2007, volume 4622 of Lecture Notes in Computer Science, pages 170-186 (Springer Berlin Heidelberg, 2007); Ravi Kannan, “Improved algorithms for integer programming and related lattice problems”, In Proceedings of the Fifteenth Annual ACM Symposium on Theory of Computing, STOC '83, pages 193-206 (New York, N.Y., USA, 1983, ACM); and Phong Q Nguyen and Damien Stehle, “Low-dimensional lattice basis reduction revisited,” 3076:338-357 (2004). A simpler version of the sampling procedure is also used in particular example implementations. The simpler version does not ensure that the distribution of procedure outcomes is close to uniform, but works well in practice. 4.4. Multiplicative Approximation Using Unit Group. In this section, the following proposition is proven: Proposition 4.6. There exists real numbers δ₀ and δ₁, . . . , δ_(d) (computed by the offline part of the procedure UNIT-ADJUST, FIG. 9) such that there exists an algorithm (online part of the procedure UNIT-ADJUST, FIG. 9) that for any real numbers t₁, . . . , t_(d) finds a unit u of

_(F) such that the following inequalities hold |log|σ_(k)(u)|−t _(k)|≤log δ_(k) for k=1, . . . ,d, under the assumption that |t₁+ . . . +t_(d)|<log δ₀. The runtime of the algorithm is bounded by a polynomial in ∥t∥.

The offline part of procedure UNIT-ADJUST computes a system of fundamental units u₁, . . . , u_(d−1)∈

_(F) ^(x) and outputs

$\delta_{k} = {\delta_{0}\sqrt{\prod\limits_{j = 1}^{d - 1}\;{\max\left\{ {{{\sigma_{k}\left( u_{j} \right)}},{{\sigma_{k}\left( u_{j}^{- 1} \right)}}} \right\}}}}$ for k=1, . . . , d, where δ₀>1 is some fixed constant.

When called with a target vector t∈

^(d) satisfying |(t, 1_(d))<log δ₀, the online part of UNIT-ADJUST simply rounds off t in the basis

$B = \begin{pmatrix} {\log{{\sigma_{1}\left( u_{1} \right)}}} & \ldots & {\log{{\sigma_{1}\left( u_{d - 1} \right)}}} & {2\;\log\;\delta_{0}} \\ \vdots & \; & \vdots & \vdots \\ {\log{{\sigma_{d}\left( u_{1} \right)}}} & \ldots & {\log{{\sigma_{d}\left( u_{d - 1} \right)}}} & {2\;\log\;\delta_{0}} \end{pmatrix}$ to the lattice vector Bm, where m=┌B⁻¹t┘∈

^(d). Then it returns the unit u=u₁ ^(m) ¹ . . . u_(d−1) ^(m) ^(d−1) . Proof that UNIT-ADJUST is Correct. Because Bm is the unique lattice vector contained in the shifted parallelepiped t+C(B), the following inequalities hold for k=1, . . . , d:

$\begin{matrix} {{{{({Bm})_{k} - t_{k}}} \leq {\max\left\{ {{x_{k}\text{:}\mspace{14mu} x} \in {\mathcal{C}(B)}} \right\}}} = {\frac{1}{2}{\max\limits_{y \in {\{{\pm 1}\}}^{d}}({By})_{k}}}} \\ {= {{\log\;\delta_{0}} + {\frac{1}{2}{\sum\limits_{j = 1}^{d - 1}\;{{\log\left( {{\sigma_{k}\left( u_{j} \right)}} \right)}}}}}} \\ {{= {\log\;\delta_{k}}},} \end{matrix}$

It is also worth noting the above shows that C(B)⊂[−log δ₁,log δ₁]× . . . ×[−log δ_(d),log δ_(d)] and ∥Bm−t∥_(B)≤1, where ∥x∥_(B):=inf{y>0:x∈C(B)y}.

Now it is shown that the running time is a polynomial in ∥t∥ and in the number of bits used to specify t_(k)

Proof that UNIT-ADJUST Runs in Polynomial Time.

Suppose that the t_(k) are given with n bits of precision. Then they can be specified using O(n+log|t_(k)|) bits as t_(k)=±2^(l−n)s, where l=┌log₂|t_(k)|┐ and s∈{0, . . . , 2^(n-1)} is an n-bit integer. First, observe that because the number field is fixed and δ₀>1 is an arbitrary fixed constant, the inverse B⁻¹ can be precomputed to sufficiently high precision and stored during the offline part. The vector m can therefore be computed in polynomial time. Also note that its norm is bounded by a polynomial in ∥t∥. Indeed,

${m} \leq {{{m - {B^{- 1}t}}} + {{B^{- 1}t}}} \leq {\frac{\sqrt{d}}{2} + {{\lambda_{\min}(B)}^{- 1}{t}}} \leq {{O\left( {t} \right)}.}$ This further implies that each |m_(i)|≤O(∥t∥), so that the output unit u=u₁ ^(m) ¹ . . . u_(d−1) ^(m) ^(d−1) can be computed by polynomially-many multiplications of the fundamental units u_(i). Therefore, log∥u∥=O(|m ₁|log∥u ₁ ∥+ . . . +|m _(d)|log∥u _(d)∥)=O(poly(∥t∥,∥u ₁ ∥, . . . ,∥u _(d)∥)), implying that the output unit can indeed be computed in polynomial time. 4.4.1. Implementation Aspects. For proving that the algorithm runs in polynomial time, it is sufficient to show that δ_(k) are fixed numbers for a given quaternionic gate set specification. It does not in principal matter how big they are. However, one can see that the additive constant C_(min) in Theorem 4.1 depends on values of δ_(k). When implementing examples of the algorithm in practice, the aim can be to achieve the smallest possible (or favorably small) constant C_(min) while maintaining good performance. For this reason, the Nearest Plane Algorithm is used instead of the simple round off procedure shown in FIG. 9. The Nearest Plane Algorithm has runtime that is polynomial in dimension of the lattice and bit sizes of the entries of vectors involved in the computation. In Section 8, the pseudo-code for an example variant of the UNIT-ADJUST procedure that has been used in practice is shown. The precision and complexity analysis above can be extended to this more practical approach. For a further analysis of the Nearest Plane Algorithm, see Phong Q. Nguyen and Damien Stehle, “An LLL algorithm with quadratic complexity,” SIAM Journal on Computing, 39(3):874-903 (2009).

The results of applying the Nearest Plane Algorithm depends on the quality of the basis used with it. In practice, and in accordance with certain example implementations, a Hermite-Korkine-Zolotarev or LLL reduction is applied to the unit lattice basis during the offline step of the algorithm. See Guillaume Hanrot and Damien Stehle, “Improved analysis of kannan's shortest lattice vector algorithm” in Alfred Menezes, editor, Advances in Cryptology—CRYPTO 2007, volume 4622 of Lecture Notes in Computer Science, pages 170-186 (Springer Berlin Heidelberg, 2007); Ravi Kannan, “Improved algorithms for integer programming and related lattice problems”, In Proceedings of the Fifteenth Annual ACM Symposium on Theory of Computing, STOC '83, pages 193-206 (New York, N.Y., USA, 1983, ACM); and Phong Q Nguyen and Damien Stehle, “Low-dimensional lattice basis reduction revisited,” 3076:338-357 (2004); Phong Q. Nguyen and Damien Stehle, “An LLL algorithm with quadratic complexity,” SIAM Journal on Computing, 39(3):874-903 (2009). This allows the procedure to further lower the contribution from log 51, . . . , log δ_(k) to the additive constant C_(min). Value of log 60 can be chosen to be very small, and its contribution to C_(min) can be made negligible without high computational overhead. Values of δ_(k) computed based on a reduced basis can also be related to the value of the regulator of number field F and known techniques for bounding the regulator can be applied to bound them.

Computing the system of fundamental units of the number field is known to be a hard problem and can be too costly even for the offline part of the algorithm. In practice, this issue can be circumvented to some extent. For example, it is sufficient to know the generators of the finite index subgroup of the unit group, but not the unit group itself. Frequently generators of such a subgroup can be computed much faster than the system of fundamental units (see, e.g., Michael E. Pohst and Hans Zassenhaus, “Algorithmic Algebraic Number Theory,” volume 30 of Encyclopedia of Mathematics and its Applications (Cambridge University Press, 1989)) or are even known in analytic form (see, e.g., Lawrence C. Washington, Introduction to Cyclotomic Fields (Springer, 1982)).

4.5. On the Conjecture Related to the Approximation Algorithm Performance.

On a high level, the performance of the example approximation algorithm depends on the properties of the set of all possible solutions to QAP (Problem 2.4). Recall the statement of QAP:

Problem 2.4.

(Quaternion approximation problem, QAP). Given

(1) A quaternion gate set specification

F, σ, a, b,

, S={

₁, . . . ,

_(M)}),

(2) target angle φ

(3) target precision ε

(4) target cost vector (L₁, . . . , L_(M)) satisfying L ₁ log N(

₁)+ . . . +L _(M) log N(

_(M))−4 log(1/ε)∈[C _(min) ,C _(max)],

where constants C_(min),C_(max) depend only on the quaternion gate set specification.

Find q from the generalized Lipschitz order (see the “generalized Lipschitz order” section) in quaternion algebra

$\left( \frac{a,b}{F} \right)$ such that nrd(q)

_(F)=

₁ ^(L) ¹ . . .

_(M) ^(L) ^(M) and ∥R_(z)(φ)−U_(q)∥≤ε, where ∥U∥=½√{square root over (Tr(UU^(†)))}.

Recall also, that map U_(q) is constructed in the “Using quaternions to represent unitaries” Section (Section 3.1) using the embedding σ:F→

that is a part of quaternion gate set specification.

Next, a formal description to the set of all solution to QAP is constructed. Let L be a generalized Lipschitz order in

$\left( \frac{a,b}{F} \right)$ (see Section 3.4). Here, the following set is used as a part of the description of all possible solutions to QAP: Sln _(r,φ,ε) ={q∈L:∥U _(q) −R _(z)(φ)∥≤ε,nrd(q)=r}

The set of all possible norms of quaternions with given cost vector L₁, . . . , L_(M) is given by nrd _(L) ₁ _(, . . . ,L) _(M) ={r∈

_(F) :r

_(F)=

₁ ^(L) ¹ . . .

_(M) ^(L) ^(M) ,r is totally positive} Using the notation above, the set of all solutions to given instance of QAP is given by:

${Sln}_{QAP} = {\bigcup\limits_{r \in {nrd}_{L_{1},\ldots\mspace{14mu},L_{M}}}{{Sln}_{r,\varphi,ɛ}.}}$ Note that for any unit u for

_(F) ^(x) it is the case that U_(q)=U_(qu). For this reason, the set of all solutions can be obtained as:

${{Sln}_{QAP} = {\bigcup\limits_{r \in {{nrd}_{L_{1},\ldots\mspace{14mu},L_{M}}/{({\mathbb{Z}}_{F}^{\times})}^{2}}}{Sln}_{r,\varphi,ɛ}}},$ where the set nrd_(L) ₁ _(, . . . , L) _(M) /(

_(F) ^(x))² is finite. Its size is equal to the size of the set: {u∈

_(F) ^(x)):u is totally positive}/(

_(F) ^(x))². Next, a discussion of the structure of the set Sln_(r,φ,ε) is provided. Consider q from Sln_(r,φ,ε). Quaternion q can be described by two elements z₁,z₂ of

_(K) as q=e ₁ ⁻¹(z ₁)+e ₂ ⁻¹(z ₂) Note that equality r=z₁z₁*−bz₂z₂ and condition ∥U_(q)−R_(z)(φ)∥≤ε imply that Re((σ_(1,+)(z ₁)−z ₀)e ^(iφ/2))≥0,σ_(k,+)(z ₁)≤σ_(k)(r) where z₀=√{square root over (σ₁(r))}(1−ε²)e^(−iφ/2). In other words z₁ belongs to the set Cand_(r,φ,ε) defined as: Cand _(r,φ,ε) ={z∈

_(K):Re((σ_(1,+)(z ₁)−z ₀)e ^(iφ/2))≥0,σ_(k,+)(z ₁)≤σ_(k)(r)}={z:σ(z)∈S _(r,φ,ε)} The observations above allows one to rewrite the set Sln_(r,φ,ε) as following:

${Sln}_{r,\varphi,ɛ}\bigcup\limits_{z_{1} \in {C\mspace{14mu}{{and}\mspace{14mu}}_{r,\varphi,ɛ}}}\left\{ {{{{e_{1}^{- 1}\left( z_{1} \right)} + {{e_{2}^{- 1}\left( z_{2} \right)}\text{:}\mspace{14mu} z_{2}}} \in {\mathbb{Z}}_{K}},{{z_{2}}^{2} = {\left( {r - {z_{1}}^{2}} \right)/\left( {- b} \right)}}} \right\}$ Note that some sets in the above union can be empty, because the relative norm equation |z₂|²=(r−|z₁|²)/(−b) does not always have a solution. Motivated by this fact, one can define the set Term_(r,φ,ε) ={z ₁ ∈Cand _(r,φ,ε): there exists z ₂∈

_(K) such that |z ₂|²=(r−|z ₁|²)/(−b)} Now assume the existence of an oracle for solving the relative norm equations and drawing points from Term_(r,φ,ε). Under this assumptions, one could have the following algorithm for solving QAP:

(1) Pick random r from nrd_(L) ₁ _(, . . . , L) _(M)

(2) Pick z₁ from Term_(r,φ,ε)

(3) Find z₂ by solving relative norm equation |z₂|²=(r−|z₁|²)/(−b)

(4) Return q=e₁ ⁻¹(z₁)+e₂ ⁻¹(z₂)

Suppose now, that such an oracle does not exist for drawing points from Term_(r,φ,ε). The algorithm can then be modified as follows:

(1) Pick random r from nrd_(L) ₁ _(, . . . , L) _(M)

(2) Pick random element z₁ from Cand_(r,φ,ε)

(3) Check if z₁ is in Term_(r,φ,ε). If this is not the case return to Step 2.

(4) Find z₂ by solving relative norm equation |z₂|²=(r−∥z₁|²)/(−b)

(5) Return q=e₁ ⁻¹(z₁)+e₂ ⁻¹(z₂)

Note that if the ratio |Term_(r,φ,ε)|/|Cand_(r,φ,ε)| were in Ω(1/log(1/ε)), then the algorithm would still run in polynomial time. In practice, there is no oracle that solves all relative norm equations in polynomial time (or even checks if there is a solution to given relative norm equation). However, if one restricts the possible right hand sides of the relative norm equation, one can check the existence of the solution and find one in polynomial time. This motivates the following definition: PolyTerm_(r,φ,ε) ={z ₁∈Term_(r,φ,ε):

((r−|z ₁|²)/(−b)) is a rational prime} This gives the example following algorithm, which is very close to the procedure APPROXIMATE in FIG. 1:

(1) Pick random r from nrd_(L) ₁ _(, . . . , L) _(M)

(2) Pick random element z₁ from Cand_(r,φ,ε)

(3) Check if z₁ is in PolyTerm_(r,φ,ε). If this is not the case return to Step 2.

(4) Find z₂ by solving relative norm equation |z₂|₂=(r−|z₁|²)/(−b)

(5) Return q=e₁ ⁻¹(z₁)+e₂ ⁻¹(z₂)

If the ratio |PolyTerm_(r,φ,ε)|/|Cand_(r,φ,ε)| were in Ω(1/log(1/ε)) and one were drawing samples from Cand_(r,φ,ε) sufficiently uniformly, the algorithm above would still run in polynomial time. In this case in the absence of the oracle for solving arbitrary norm equation. Note the above discussion implies, that

$\frac{{PolyTerm}_{r,\varphi,ɛ}}{{Cand}_{r,\varphi,ɛ}} = {{\frac{{PolyTerm}_{{ru}^{2},\varphi,ɛ}}{{Cand}_{{ru}^{2},\varphi,ɛ}}\mspace{14mu}{for}\mspace{14mu}{any}\mspace{14mu} u} \in {{\mathbb{Z}}_{F}^{\times}.}}$ For this reason, the ratio above is well defined for r/(

_(F) ^(x))². The conjecture that implies that the algorithm runs in polynomial time is the following: Conjecture 4.8. Keeping the notation introduced before in this section, for any r from nrd _(L) ₁ _(, . . . ,L) _(M) /(

_(F) ^(x))² ={r∈

_(F) :r

_(F)=

₁ ^(L) ¹ , . . . ,

_(M) ^(L) ^(M) ,r is totally positive}/(

_(F) ^(x))² the ratio |PolyTerm_(r,φ,ε)|/|Cand_(r,φ,ε)| is in Ω(1/log(1/ε)).

5. Relative Norm Equations

In this section, it is shown how a solution z∈K to a relative norm equation of the form N_(K/F)(z)=e between a CM field K and its totally real subfield F=K∩

can be efficiently computed, provided such solutions exist at all for the given right hand side e∈

_(F). The totally positive element e arises from the RANDOM-INTEGER-POINT step in lines 9 and 10 in the main algorithm (see FIG. 1). Formally, this concerns the following problem:

Problem 5.1

(CM Relative Norm Equation). Let K/F be a CM field of constant degree over

and let e be a totally positive element of

_(F). The task is to find an element z of

_(K) such that N_(K/F)(z)=zz*=e in time polynomial in the bit-size of e, provided such an element z exists.

In the following paragraph, an example approach to solving a relative norm equation as in Problem 5.1 is described and pseudo-code implementations of IS-EASILY-SOLVABLE step in line 11 and FAST-SOLVE-NORM-EQ in line 12 of the main algorithm in FIG. 1 are given.

Relative norm equations of the form N_(K/F)(z)=e have been studied in the literature before. Early approaches include various methods that proceed by establishing a bounding box that will contain a solution provided it exists and then checks the candidates in the bounding box. See, e.g., Zenon I. Borevich and Igor R. Shafarevich, “Number theory,” volume 20 of Pure and Applied Mathematics. Academic Press (1967); Claus Fieker, Andreas Jurk, and Michael E. Pohst, “On solving relative norm equations in algebraic number fields,” Math. Comput., 66(217):399-410 (1997); and Dennis A. Garbanati, “An algorithm for finding an algebraic number whose norm is a given rational number,” J. Reine Angew. Math., 316:1-13, 1980; Michael E. Pohst and Hans Zassenhaus, Algorithmic Algebraic Number Theory, volume 30 of Encyclopedia of Mathematics and its Applications (Cambridge University Press, 1989). Unfortunately, these methods are exponential in the bit-size of the right hand side. Next, there is a method based on S-units. See, e.g., Henri Cohen, “Advanced Topics in Computational Number Theory,” Graduate Texts in Mathematics (Springer, New York, 2000); and Denis Simon, “Norm equations in relative number fields using S-units,” Mathematics of Computation, 71(239):1287-1305 (2002). This requires the factorization of the right hand side of the equation, along with precomputation of the relative class group of the extension K/F, and some additional data that is dependent on the right hand side. Therefore, it is not clear that the resulting algorithm runs in polynomial time. Relative norm equations have also been studied in the context of cryptanalysis of lattice based cryptography, e.g., of the NTRU system. Further, in Craig Gentry and Michael Szydlo, “Cryptanalysis of the Revised NTRU Signature Scheme,” in International Conference on the Theory and Applications of Cryptographic Techniques (EUROCRYPT'02), pages 299-320 (2002), an algorithm is described to solve relative norm equations for cyclotomic fields over their totally real subfield. This algorithm uses Fermat's little theorem for ideals in

_(K) in conjunction with LLL reduction to find a solution, which is known to exist in the context in which the algorithm is applied. See also Sanjam Garg, Craig Gentry, and Shai Halevi, “Candidate multilinear maps from ideal lattices,” in International Conference on the Theory and Applications of Cryptographic Techniques (EUROCRYPT'13), pages 1-17 (2013). However, like Simon's S-unit based algorithm, the algorithm relies on some properties of the right hand side, and therefore does not seem to run in time that is a polynomial in the bit-size of the right hand side.

In embodiments of the disclosed technology, a different route is taken in which precompute a finite set of attributes of K/F that does not depend on the right hand side e. The example method reduces the problem for general right hand side to a bounded size instance. Furthermore, the approach leverages the fact that K is a CM field and that the right hand side is of a particular form, which are called benign integers. These are characterized in terms of a finite set

={

₁, . . . ,

_(k)}⊂

_(F):

Definition 5.2

(

-Benign Integers). Let

be a set of prime ideals of

_(F). An integer e∈

_(F) is called benign if it is totally positive and the prime factorization of the ideal generated by e satisfies

$\begin{matrix} {{{e\;{\mathbb{Z}}_{F}} = {{\mathfrak{q}}{\prod\limits_{{\mathfrak{p}} \in \mathcal{P}}\;{\mathfrak{p}}^{e{({\mathfrak{p}})}}}}},} & (16) \end{matrix}$ where q is prime and e(

)≥0 for all

∈

.

The primes in

are defined by the user parameter P in algorithm APPROXIMATE (FIG. 1) and are precomputed by the offline part of procedure IS-EASILY-SOLVABLE (FIG. 4).

5.1. Measuring the Bit-Size of the Input.

There are several natural ways to measure the bit-size of the algebraic numbers that are involved as the input and the output of a relative norm equation. These definitions are briefly discussed and it is shown that in the example approach they are all within a constant factor of each other.

Let K/

be a Galois extension and let B={b₁, . . . , b_(n)} be a basis for

_(K) over

, where n=[K:

], i.e.,

_(K)=⊕_(i=1) ^(n)b_(i)

. Any x∈

_(K) can then be represented as x=Σ_(i=1) ^(n)b_(i)x_(i) where x_(i)∈

, e.g., one can define the bit-size with respect to B as ∥x∥_(B):=Σ_(i=1) ^(n)|x_(i)|. Alternatively, for CM fields, have n=2d where [F:

]=d and one can use the quadratic form T₂(x)=Σ_(i=1) ^(d)|σ_(i)(x)|² as a measure for the bit-size of x. Also, one can use a notion of bit-size that is valid for general ideals I∈

_(K) and not just for the principal ideals e

_(F): following Jean-Francois Biasse and Claus Fieker, “A polynomial time algorithm for computing the HNF of a module over the integers of a number field,” in International Symposium on Symbolic and Algebraic Computation (ISSAC'12), pages 75-82 (2012), a matrix M∈

^(n×n) is chosen for a basis of I expressed over an integral basis of

_(K). If M is in Hermite Normal Form, then each entry can be bounded by |det(M)|=

(I), e.g., one can define S(I):=n² log₂(

(I)) as the bit-size of I. For principal ideals, the following definition is applied S(x):=n log₂(max_(i)|x_(i)|).

It turns out that T₂(x) and S(x) are related. More precisely, the following result can be found (see, e.g., Jean-Francois Biasse and Claus Fieker, “A polynomial time algorithm for computing the HNF of a module over the integers of a number

field,” In International Symposium on Symbolic and Algebraic Computation (ISSAC'12), pages 75-82 (2012)):

Lemma 5.3.

Let K/

be a CM field, let x∈

_(K), and let Δ_(K)=det(T₂(b_(i),b_(j)))² denote the discriminant of K. Then the bound ½ log₂(T₂(x))≤Õ(S(x)/n+n²+log₂(Δ_(K))) holds. Furthermore, one has that S(x)≤Õ(d(d+½ log₂(T₂(x)))).

Next, a bound is established that allows to relate T₂(x) to the bit-size of the expansion of x with respect to any given basis B.

Lemma 5.4.

Let K/

be a CM field and let B={b₁, . . . , b_(n)} be an integral basis for K over

, where n=[K:

]. For x=Σ_(i=1)x_(i)b_(i)∈K define the bit size of x with respect to the basis B as ∥x∥_(B):=Σ_(i=1) ^(n)|x_(i)|². Let M=[Tr(b_(i)b_(j)*)]_(i,j=1, . . . , n) be the Gram matrix of B and let λ_(max) and λ_(min) be the largest, respectively smallest eigenvalue of M. Then λ_(min) ∥x∥ _(B)/2≤T ₂(x)≤λ_(max) ∥x∥ _(B)/2. Proof. Let x∈K and let x=Σ_(i=1) ^(n)x_(i)b_(i) be its expansion over the chosen basis. Recalling that T₂(x)=

(xx*)=

(xx*)/2 one can obtain that T₂(x)=Σ_(i,j=1) ^(n)x_(i)x_(j)*Tr(b_(i)b_(j)*)/2. One can rewrite this as T₂(x)=(x₁, . . . , x_(n))M(x₁, . . . , x)^(t)/2 where M is the integer valued, symmetric, and positive-definite matrix with entries M_(i,j)=Tr(b_(i)b_(j)*). By diagonalizing M in an eigenbasis, one sees that λ_(min)∥x∥_(B)/2≤T₂(x)≤λ_(max)∥x∥_(B)/2 as claimed. □

In an example approach taken in this disclosure, the field K is considered to be a constant. This implies that quantities such as the degree [K:

] or the discriminant Δ_(F) of the totally real subfield F=K⁺ are constants. By choosing B to be an LLL-reduced basis, one can obtain from Lemma 5.4 that T₂(x) and ∥x∥_(B) are related by a constant factor. See also Karim Belabas, “Topics in computational algebraic number theory,” Journal de Theorie des Nombres, 16:19-63 (2004). To summarize, Lemmas 5.3 and 5.4 imply that all measures of bit-size considered in the following are within constant factors of each other.

5.2. Hermite Normal Form and Lifting of Ideals.

In order to represent ideals in rings of integers, Hermite normal forms are a useful tool. For example, see the result from Arne Storjohann and George Labahn, “Asymptotically fast computation of Hermite normal forms of integer matrices,” in Proceedings of the 1996 International Symposium on Symbolic and Algebraic Computation, ISSAC '96, pages 259-266 (1996) that allows to give a polynomial bound on both, the complexity of computing a Hermite Normal Form (HNF) of an integer matrix, and the bit-size of the output. See also Arne Storjohann, “Algorithms for matrix canonical forms,” PhD thesis, ETH Zurich (2000) for a discussion and comparison with other efficient algorithms to compute HNFs. Theorem 5.5. Let A∈

^(n×m) be a rank r integer matrix and let ∥A∥:=max_(i,j)|A_(i,j)|. There exists a deterministic algorithm that computes the HNF of A in time Õ(m^(θ)n log∥A∥), where the Õ notation ignores log-factors and 2≤θ≤2.373 is the exponent for matrix multiplication.

Next, a discussion will be presented of how to lift primes ideals in

_(F) to prime ideals in

_(K). Write K=F(θ), where θ is a primitive element. Recall that the conductor of K/F is defined as

={x∈

_(K):x

_(K)⊆

_(F)[θ]}. It is well known that at least the prime ideals that are coprime with

can be easily lifted via a reduction to factoring the minimal polynomial of θ over a suitable finite field as described below. See, e.g., Jurgen Neukirch, “Algebraic Number Theory”, Ch. I, Prop. 8.3 (Springer, 1999); Daniel A. Marcus, “Number Theory,” Ch. 3, Thm. 27 (Springer 1977); or Henri Cohen, “Advanced Topics in Computational Number Theory,” Graduate Texts in Mathematics, Prop. 2.3.9 (Springer New York, 2000).

Being coprime with the conductor is equivalent to being coprime with

_(K)/

_(F)[θ]| which leaves only a finite set of primes for which it does not hold, e.g., these primes will be added to the exceptional set

.

Theorem 5.6.

Let K/F be a Galois extension where K=F(θ). Denote by g(X)∈F[X] the minimal polynomial of θ and by

the conductor of K/F. Let

be a prime ideal in

_(F) that is coprime with

. Let

=

_(F)/

be the finite field corresponding to the residues mod

and let g (X)= g ₁(X)^(e) ¹ . . . g _(r)(X)^(e) ^(r) ∈

[X] be the factorization of g(X)=g(X) mod

into irreducible polynomials g _(i)(X)=g_(i)(X) mod

∈

[X] where all g _(i)(X) are chosen to be monic polynomials. Then the ideals

_(i):=

_(K) +g _(i)(θ)

_(K)  (17) are precisely the prime ideals of

_(K) that are lying over

and all these ideals are pairwise different.

Next, the time-complexity of computing the list of ideals

_(i) lying over

is analyzed. Factoring of a polynomial f(X) of degree n over a finite field

_(q) of size q is known to run in time polynomial in n and log(q). See, e.g., Joachim von zur Gathen and Jurgen Gerhard, “Modern Computer Algebra,” Theorem 14.14 (Cambridge University Press, 2nd edition, 2003).

Theorem 5.7.

Let

_(q) be a finite field and let f(X)∈

_(q)[X] be a polynomial of degree n. Then there exists a probabilistic polynomial time algorithm that computes the factorization of f(X)=Π_(i=1) ^(k)f_(i)(X) into irreducible polynomials over

_(q). The probability of success of the algorithm is at least ½ and the expected running time can be bounded by Õ(n² log q).

The algorithm in Theorem 5.7 proceeds in 3 stages, namely (i) squarefree factorization, (ii) equal degree factorization, and (iii) distinct degree factorization. As in the present case [K:F]=2, there are only two possibilities of possible splitting behavior of g(X): either a) this polynomial is irreducible which according to Theorem 5.6 corresponds to the case in which

is inert as g_(i)(θ)=0 and

itself generates a prime ideal in

_(K) or b) the polynomial splits as g(X)=g ₁(X)g ₂(X) where both factors are linear. This corresponds to the case where there are two ideals

₁ and

₂ lying over

. On account of Galois theory in this case,

₂=

₁*.

Using Theorem 5.7, one can obtain a refined version of Theorem 5.6 that bounds the running-time of finding the ideals lying over a given prime ideal

in terms of the bit complexity of

.

Corollary 5.8.

Let K/F be a CM field where K=F(θ). Let

⊆

_(F) be a prime ideal and let n=S(

) be the bit-size of

. Then there exists a polynomial time algorithm to compute all ideals

lying over

. Furthermore, the bit-size of the

is polynomial in n.

Proof.

Let g(X)∈F[X] be the minimal polynomial of θ and let

_(F):=O_(F)/

be the finite field that arises as the residue field of

. Using Theorem 5.7, one sees that the factorization g(X)=Π_(i=1) ^(r) g _(i)(X) can be computed into irreducibles in polynomial time in n. For each factor g _(i), a corresponding lifting g_(i)(X) can be determined by considering the components in of

_(q) as elements of F while maintaining the same bit-size, e.g., one can obtain that the element g_(i)(θ) has a bit-size that is polynomial in n as well. Now, the 2-generator representation of

is generated as in Eq. (17). Next, an HNF representation of this ideal is found: as θ obeys a quadratic equation F one can choose B={1, θ} as a basis of K/F. Hence, one can obtain a set of generators in the form

$H = {\begin{pmatrix} A & 0 & a & b \\ 0 & A & c & d \end{pmatrix} \in {\mathbb{Z}}^{2\; d \times 4\; d}}$ where n=2d, A is an HNF for

, and the matrix

$\quad\begin{pmatrix} a & b \\ c & d \end{pmatrix}$ is the expansion over

of the linear map that describes the multiplication by g_(i)(θ) with respect to B. Note that all coefficients of H are bounded in bit-size by a polynomial in n. Theorem 5.5 can be used to compute an HNF H′ for H, and hence the ideal

, in time polynomial in n. Theorem 5.5 also implies that the output size, e.g., all coefficients of H′, are polynomial in n. 5.3. Outline of an Example Algorithm to Solve Relative Norm Equations. Let K be a CM field over its totally real subfield F:=K⁺ and denote by {σ₁, . . . , σ_(d)} the real embeddings of F into

, e.g., [K:F]=2 and [F:

]=d. The extension K/F is Galois and its Galois group is generated by complex conjugation, e.g., Gal(K/F)=

τ

, where τ(x):=x*. Furthermore, denote by

_(F) and

_(K) the rings of integers in K and F, respectively. Recall that for an ideal I⊆

_(K) the norm is defined as N_(K/F)(I):=|

_(K)/I| which for principal ideals I=(x) coincides with the usual definition as the product of all Galois conjugates, i.e., N_(K/F)((x))=x·x*.

For a given element e ∈

_(F), one example approach to solving the relative norm equation N_(K/F)(z)=e where z∈

_(K) relies on the observation that if I·I*=e

_(F) is a factorization of ideals in

_(F) and η is an arbitrary non-zero element in the lattice generated by I, then N_(K/F)(I)|N_(K/F)(η). This alone would not be a very useful property as potentially the quantity on the right might be unbounded. Using the fact that [K:F]=2 is constant and that K is a CM field, however, one will be able to show that for suitable η the quotient N_(K/F)(η)/N_(K/F)(I) will be a constant that just depends on the CM field K and not on the right hand side e of the norm equation (5.1). Furthermore, it will be shown that one can find such z in polynomial time, provided that e is benign.

In the next two sections, it is shown that there exists a probabilistic algorithm that runs in polynomial time (with respect to the bit-size of the right hand side e) and finds an element z such that N_(K/F)(z)=e or else reports that no such element exists. This is first shown for the case where e generates a prime ideal and then in a subsection section for the case of general benign e. The example algorithm proceeds in several stages:

-   Step 1 Determine whether the right hand side e     _(F) is prime, respectively benign. -   Step 2 If e     _(F) is prime, respectively benign, extract     from the prime decomposition eq. (16). -   Step 3 Compute a prime ideal     lying over     using Theorem 5.6. -   Step 4 For all ideals     _(i) in     , precompute all ideals     _(i,j) lying over them. This can be done offline. -   Step 5 For all exponents e_(i,j), where the tuple e is taken from a     bounded set of candidates, compute I:=     ⋅     _(1,1) ^(e) ^(1,1) ⋅ ⋅ ⋅ . . . ⋅     _(k)     and apply the following steps to all candidates. In the case of     prime e     _(F) this step can be omitted and I:=     . -   Step 6 Compute LLL reduced basis for lattice I to obtain an     approximation η to the SVP for I. -   Step 7 Compute the norm of η and compute γ:=N_(K/F)(η)/N_(K/F)(I)     and attempt to solve the norm equation for γ using a known method,     such as the method described in Denis Simon, “Norm equations in     relative number fields using S-units,” Mathematics of Computation,     71(239):1287-1305 (2002). Let w be such that N_(K)/F(w)=γ. -   Step 8 Output “fail” if no such solution exists. Otherwise, return a     solution z=ηw∈     _(K).

Provided below is a proof that all steps can be performed by a classical algorithm whose runtime is polynomial in the bit-size of e. In the “Constructing the solution: prime case” Section (Section 5.5), this is shown for the somewhat simpler, however in practice frequently occurring, case where the right hand side e

_(F) generates a prime ideal. Further, in the “Constructing the solution: general case” Section (Section 5.6), a discussion is provided of how the case of any benign e can be handled. Before proceeding with the proof, however, it is helpful to show another technical result, namely that it is indeed possible to find an element r as needed for Steps 6 and 7 such that the co-factor γ:=N_(K/F)(η)/N_(K/F)(I) is bounded.

5.4. Bounding the Co-Factor.

Assuming that the norm equation I·I*=e

_(K) is solvable implies that I=ξ

_(K) for some element ξ∈

_(K). It is now shown that one can find an element η∈e

_(K) such that the quotient of the norms of I and η

_(K) is a constant. Here, I is considered to be a 2d-dimensional

lattice, where d=[F:

]. This means that there exists a basis {a₁, . . . , a_(2d)}⊂K such that I=a₁

+ . . . +a_(2d)

. Recall further that there is a quadratic form on I defined by (x,y):=

(xy*) and that the Gram matrix G_(i,j):=(a_(i),a_(j)) is integer valued, i.e., G∈

^(2d×2d). Furthermore, for the volume of the fundamental parallelepiped of I, the identity vol(I)=√{square root over (det(G))} holds. See, e.g., Daniele Micciancio and Sha Goldwasser, “Complexity of lattice problems: a cryptographic perspective,” volume 671 of Springer International Series in Engineering and Computer Science (Springer, 2002).

The fact that I=ξ

_(K) is principal is used in the following lemma to compute vol(I) in terms of the absolute norm of ξ:

Lemma 5.9.

For each I=ξ

_(K), vol(I)=vol(

_(K))·

(ξ).

Proof.

First, a basis {b₁, . . . , b_(2d)} of

_(k) over

, i.e.,

_(K)=b₁

+ . . . +b_(2d)

is chosen. With respect to this basis, multiplication with the fixed element ξ is a linear transformation M_(ξ) defined via M_(ξ)(x₁, . . . , x_(2d))=ξ(b₁x₁+ . . . +b_(2d)x_(2d)), and the determinant of M_(ξ) is equal to the norm

(ξ). Note also, as K is a CM field, all Galois automorphisms σ∈Gal(K/

) come in complex conjugate pairs, i.e.,

(ξ)≥0, i.e.,

(ξ)=|det(M_(ξ))|.

By applying a base change to the Gram matrix G in which one goes from pairs of conjugates σ_(i),σ _(i) to Re(σ_(i)), Im(σ_(i)), where by convention the Galois automorphisms are ordered in such a way that the first d are pair-wise non-conjugates under complex conjugations, e.g., σ _(i)≠σ_(j) for all 1≤i<j≤d.

For general x_(i)∈K one can denote by V(x₁, . . . , z_(2d)) the matrix

${V\left( {x_{1},\ldots\mspace{14mu},x_{2\; d}} \right)} = {\begin{pmatrix} {{Re}\left( {\sigma_{1}\left( x_{1} \right)} \right)} & \ldots & {{Re}\left( {\sigma_{1}\left( x_{2\; d} \right)} \right)} \\ {{Im}\left( {\sigma_{1}\left( x_{1} \right)} \right)} & \ldots & {{Im}\left( {\sigma_{1}\left( x_{2\; d} \right)} \right)} \\ \vdots & \ddots & \vdots \\ {{Re}\left( {\sigma_{d}\left( x_{1} \right)} \right)} & \ldots & {{Re}\left( {\sigma_{d}\left( x_{2\; d} \right)} \right)} \\ {{Im}\left( {\sigma_{d}\left( x_{1} \right)} \right)} & \ldots & {{Im}\left( {\sigma_{d}\left( x_{2\; d} \right)} \right)} \end{pmatrix}.}$

Using this matrix, one can then express the volume of I as vol(I)=√{square root over (det(G))}=|det(V(a₁, . . . , a_(2d)))|, where the set {a_(i):i=1, . . . , 2d} forms a basis for I over

.

Next, one can observe that the matrices V_(I):=V(a₁, . . . , a_(2d)) and V

_(K) :=V(b₁, . . . , b_(2n)) are related via V_(I)=M_(ξ)

_(K). From this, one can conclude that vol(I)=√{square root over (det(G))}=|det(V _(I))|=|det(M _(ξ))|·|det(V

_(K) )=

(ξ)·vol(

_(K)), as claimed.

It is now shown how to find an element η∈ξ

_(K) such that the quotient of the norms of I and η

_(K) is a constant as mentioned in the beginning of this section.

Lemma 5.10.

Let I=ξ

_(K) be an ideal in

_(K) such that I·I*=e

_(K). Then there exists η∈

_(K) such that

(ξ/η) is upper bounded by a constant C_(K) that depends just on the extension K/

.

Proof.

As above, consider I as a lattice I=a₁

+ . . . +a_(2d)

. Now use the LLL algorithm on the basis {a₁, . . . , a_(2d)}⊂

^(2d). Using the LLL algorithm described in Phong Q. Nguyen and Damien Stehle, “Floating-point LLL revisited,” in 24th Annual International Conference on the Theory and Applications of Cryptographic Techniques (EUROCRYPT'05), pages 215-233 (2005) it is known that the first vector v=b₁ in the LLL reduced basis {b₁, . . . , b_(2d)} for I satisfies the following bound: ∥v∥≤(4/3)^((2d−1)/4)vol(I)^(1/2d). As I=ξ

_(K) is by assumption principal and v∈I, there exists an element η∈

_(K) such that v=ηξ. In order to finish the proof of the lemma, it remains to show that the norm of η is upper bounded by a constant that just depends on K alone and is in particular independent of the right hand side e of the norm equation:

${N_{K/{\mathbb{Q}}}(\eta)} = {{{N_{K/{\mathbb{Q}}}(v)}/{N_{K/{\mathbb{Q}}}(\xi)}} = {{{\prod\limits_{i = 1}^{d}\;{{{\sigma_{i}(v)}}^{2}/{N_{K/{\mathbb{Q}}}(\xi)}}}\overset{AGM}{\leq}{\left( {\frac{1}{d}{\sum\limits_{i = 1}^{d}\;{{\sigma_{i}(v)}}^{2}}} \right)^{d}/{N_{K/{\mathbb{Q}}}(\xi)}}} = {{\left( {\frac{1}{2\; d}{{Tr}_{K/{\mathbb{Q}}}\left( {vv}^{*} \right)}} \right)^{d}/{N_{K/{\mathbb{Q}}}(\xi)}} = {{{\left( \frac{1}{2\; d} \right)^{d}{{v}^{2\; d}/{N_{K/{\mathbb{Q}}}(\xi)}}}\overset{LLL}{\leq}{\left( {\frac{1}{2\; d}\left( \frac{4}{3} \right)^{{({{2\; d} - 1})}/2}} \right)^{d}{{{vol}(I)}/{N_{K/{\mathbb{Q}}}(\xi)}}}}\overset{{Lemma}\mspace{14mu} 5.9}{=}{{\left( {\frac{1}{2\; d}\left( \frac{4}{3} \right)^{{({{2\; d} - 1})}/2}} \right)^{d}{N_{K/{\mathbb{Q}}}(\xi)}{{{vol}\left( {\mathbb{Z}}_{K} \right)}/{N_{K/{\mathbb{Q}}}(\xi)}}} = {{\left( {\frac{1}{2\; d}\left( \frac{4}{3} \right)^{{({{2\; d} - 1})}/2}} \right)^{d}{{vol}\left( {\mathbb{Z}}_{K} \right)}} = {:C_{K}}}}}}}}$ where the first inequality is the arithmetic-geometric-mean inequality (AGM) and v is the first basis vector obtained via LLL reduction for δ_(LLL)≡1 and η_(LLL)≡½ as in Phong Q. Nguyen and Damien Stehle, “Floating-point LLL revisited,” in 24th Annual International Conference on the Theory and Applications of Cryptographic Techniques (EUROCRYPT'05), pages 215-233 (2005).

Using Lemma 5.10 will finally put one into a position to solve norm equations as in Problem 5.1 efficiently, in case e is a benign number as defined in Definition 5.2. In the “Constructing the solution: prime case” Section (Section 5.5) this result is established first for the case where e

_(F) is a prime ideal itself as this case is relatively straightforward. Then, the more general case of benign e will be discussed in the “Constructing the solution: general case” Section (Section 5.6).

5.5. Constructing the Solution: Prime Case.

A proof of the following theorem is now provided, which is restated from the “High level description of the algorithm” Section (Section 4.1).

Theorem 4.4.

Given totally positive element e of F, there exists an algorithm for testing if the instance of integral relative norm equation in K/F zz*=e,z∈

_(K) can be solved in polynomial time in log T₂(e) (procedure IS-EASILY-SOLVABLE, where T₂(e)=Σ_(i=1) ^(d)σ_(k) ²(e)). If the test is passed, there exist another algorithm for deciding if the solution exists and finding it that runs in time polynomial in log T₂(e) (procedure FAST-SOLVE-NORM-EQ). Procedure IS-EASILY-SOLVABLE returns true for at least those cases when the ideal e

_(F) is prime.

Now, with all of these pieces in place, it possible to prove that the norm equations arising in the context of the example approximation method can be solved in time that is polynomial in the input size. First, proof is provided for the case where the right hand side e generates a prime ideal

=e

_(F), then, proof is provided for the more general case of a benign integer e for the next section.

Proof (of Theorem 4.4).

Let n=S(e) be the bit-size of the ideal generated by e. Next, a discussion will be provided for all Steps 1-7 described in Section 5.3 and it will be shown that all operations can be performed in time that is upper bounded by a polynomial in n.

In Step 1, one can run a test which is described in subroutine IS-EASILY-SOLVABLE shown in FIG. 5. The norm computations in lines 6 and 7 are clearly polynomial in n. Finally, for the primality test at the last line of the subroutine one can use for instance a probabilistic test such as Miller-Rabin or a deterministic test. See, e.g., Henri Cohen, “A Course in Computational Algebraic Number Theory,” Graduate Texts in Mathematics (Springer, 1993) and Manindra Agrawal, Neeraj Kayal, and Nitin Saxena, “PRIMES,” in P. Annals of Mathematics, 160(2):781-793 (2004).

Step 2 can be done by computing quotients of the form (e)(

_(i))⁻¹ which can be done in polynomial time in the input bit-size and at an increase per division that is also at most polynomial. See, e.g., Jean-Francois Biasse and Claus Fieker, “A polynomial time algorithm for computing the HNF of a module over the integers of a number field” in International Symposium on Symbolic and Algebraic Computation (ISSAC'12), pages 75-82 (2012). Eventually this yields the prime ideal

.

This step is also done in subroutine IS-EASILY-SOLVABLE. All subsequent steps and line numbers refer to subroutine FAST-SOLVE-NORM-EQ shown in FIG. 5.

For Step 3, in line 12 of subroutine FAST-SOLVE-NORM-EQ, Theorem 5.6 can be used and the complexity analysis given in Corollary 5.8 can be used in order to compute an HNF for the ideal

lying over

in polynomial time and almost polynomial increase of the bit-size.

Steps 4 in line 3 is an offline computation which does not count toward the cost of the online solution of the norm equation.

Step 5 does not have to be carried out, as by assumption in this subsection it is assumed that

is prime (e.g., there is only one prime ideal

to consider). This step and the consequences for the subsequent steps in case e is benign but not prime are discussed in the next subsection.

Step 6 in line 13 involves the computation of a reduced lattice basis for the ideal corresponding to

from the HNF that was computed in Step 3. Using bounds on the complexity of the LLL algorithm one can see that the running time of this step is polynomial in the input size n and so is the bit-size of the short vector η that is produced by this computation. See, e.g., Phong Q. Nguyen and Damien Stehle, “An LLL algorithm with quadratic complexity,” SIAM Journal on Computing, 39(3):874-903 (2009).

For Step 7 in line 15, a method for solving norm equations can be used, such as Simon's S-unit based algorithm for which an implementation is available. See, e.g., Henri Cohen, “A Course in Computational Algebraic Number Theory,” Graduate Texts in Mathematics (Springer, 1993). As the element γ is constant and does not depend on the input size n one can assume that this computation can be done in constant time that does not affect the overall running time.

Finally, in Step 8, the case where there is no solution for γ is considered, which then implies that there is no solution for e and combine η and w using an ideal multiplication into the final solution z to the norm equation N_(K/F)(z)=e. This again can be done in polynomial time.

A pseudocode description of the Steps 3-8 is given in subroutine FAST-SOLVE-NORM-EQ shown in FIG. 5. From the above discussion of the steps it follows that the overall runtime is polynomial in the input size n=S(e) and by Lemma 5.3 therefore also polynomial in the bit-size log T₂(e) as claimed.

5.6. Constructing the Solution: General Case.

Here, a brief discussion is provided of the implications of e being benign but not prime. In certain example implementations, this will involve a change in Step 5; for example, instead of only considering the prime ideals

lying over

, the procedure considers all ideals I that can be formed by multiplying

with the ideals lying over the various prime factors

_(i)∈

. All ideals lying over

can be precomputed without any additional cost to the online part. Also, it is noted that even though the number of ideals to be considered grows significantly, this increase is still just a constant as for any given input parameter P of FAST-SOLVE-NORM-EQ in FIG. 5 the size of

=

₀∪

₁ is constant. The prime ideals in

₀ can be found by enumeration and those in

₁ can be found by factoring the discriminant Δ_(F) which, as argued before, can be

done in polynomial time as F is constant. In order to compute I, one can now multiply

with all possible combinations of factors that are extensions

of ideals in

. Each such multiplication can be done in polynomial time with polynomial increase in bit-size. See Jean-Francois Biasse and Claus Fieker, “A polynomial time algorithm for computing the HNF of a module over the integers of a number field,” in International Symposium on Symbolic and Algebraic Computation (ISSAC'12), pages 75-82 (2012). Then one can perform each of the subsequent Steps 6-7 as in the prime case and report success if any of the considered cases leads to a successful solution z and “failure” otherwise. Overall, also in case of benign right hand sides e, one can obtain a polynomial time classical algorithm to solve the norm equation N_(K/F)(z)=e. 5.7. Performance Improvement: Filtering Out Candidates Via the Hasse Principle. It is possible to perform tests as to whether a solution to N_(K/F)(z)=e over the rational elements of K (not necessarily elements in

_(K)) exists. A known test in this regard is the Hasse Norm Theorem that asserts that a global solution, e.g., a solution over K, exists if and only if a solution exists with respect to all local fields associated with K/F. See, e.g., Jurgen Neukirch, “Algebraic Number Theory,” (Springer, 1999). More precisely: Theorem 5.11 (Hasse Norm Theorem). Let K/F be a cyclic extension. An element e∈F^(x) is a norm of an element in K^(x) if and only it is a norm at every prime of F, including the infinite primes.

In practice, it is not necessary to check all primes of F, a finite set of primes is sufficient: as described in Vincenzo Acciaro and Jurgen Kluners, “Computing local Artin maps, and solvability of norm equations,” J. Symb. Comput., 30(3):239-252 (2000), the only primes that need to be checked are a) the divisors of the conductor

of K/F and b) all finite primes dividing the ideal e

_(F). If e is benign, one can therefore efficiently compute the prime factorization and hence can perform this sufficient test for solvability of the norm equation. Note that this test can only be used in this one-sided sense as there are examples of degree 2 extensions K/F known where for e∈

_(F) the equation N_(K/F)(z)=e is solvable over K but not over

_(K). In practice, the test is reasonably fast in order to eliminate some candidates e. Example pseudo-code for this test PASS-HASSE-SOLVABILITY-TEST is summarized in FIG. 6.

6. End to End Examples of Using the Framework

In this section, examples of using an exemplary implementation of the disclosed technology are provided. In two examples, an example implementation of the disclosed framework is applied to reproduce results on Clifford+T and V-basis gate sets. The exact synthesis framework disclosed in Vadym Kliuchnikov and Jon Yard, “A framework for exact synthesis,” arXiv:1504.04350 (2015); U.S. Provisional Patent Application No. 62/146,182 filed on Apr. 10, 2015, and entitled “METHOD AND SYSTEM FOR QUANTUM CIRCUIT SYNTHESIS USING QUATERNION ALGEBRA”; and PCT International Application No. PCT/US2016/025958 filed on Apr. 5, 2016, and entitled “METHOD AND SYSTEM FOR QUANTUM CIRCUIT SYNTHESIS USING QUATERNION ALGEBRA” is used for the end to end compilation. For this reason, this discussion describes how to use this framework to reproduce previously known results on exact synthesis from Alex Bocharov, Yuri Gurevich, and Krysta M. Svore, “Efficient decomposition of single-qubit gates into V basis circuits,” Physical Review A, 88(1):1-13 (July 2013); Vadym Kliuchnikov, Dmitri Maslov, and Michele Mosca, “Fast and efficient exact synthesis of single qubit unitaries generated by Clifford and T gates,” Quantum Information and Computation, 13(7-8):0607-0630 (June 2013); and Vadym Kliuchnikov, Dmitri Maslov, and Michele Mosca, “Asymptotically optimal approximation of single qubit unitaries by Clifford and T circuits using a constant number of ancillary qubits,” Physical Review Letters, 110(19):1-5 (December 2012), and also results on exact synthesis over Clifford+R_(z)(π/16) (see Simon Forest, David Gosset, Vadym Kliuchnikov, and David McKinnon, “Exact synthesis of single-qubit unitaries over clifford-cyclotomic gate sets,” Journal of Mathematical Physics, 56(8):082201 (2015)).

Two other examples correspond to approximating using gate sets Clifford+R_(z)(π/16) and Clifford+T+V. No number-theoretic style approximation algorithms for this gate sets were known before. A brief summary of the examples is given in Table 1. Amongst other data, the table contains time needed for precomputation stage for all of the disclosed examples and the value of the additive constant appearing in Theorem 4.1. All the data about the example algorithm performance is based on an example implementation of exact synthesis (see, e.g., Vadym Kliuchnikov and Jon Yard, “A framework for exact synthesis,” arXiv:1504.04350 (2015); U.S. Provisional Patent Application No. 62/146,182 filed on Apr. 10, 2015, and entitled “METHOD AND SYSTEM FOR QUANTUM CIRCUIT SYNTHESIS USING QUATERNION ALGEBRA”; and PCT International Application No. PCT/US2016/025958 filed on Apr. 5, 2016, and entitled “METHOD AND SYSTEM FOR QUANTUM CIRCUIT SYNTHESIS USING QUATERNION ALGEBRA”.) and approximation frameworks using the computer algebra system MAGMA. The total number of lines of code needed for it implementation is about 2500. More examples of running approximation stage of the algorithm are provided in Section 7.

6.1. Clifford+T.

In this section, it is described how to obtain results from Peter Selinger, “Efficient Clifford+T approximation of single-qubit operators,” (December 2012) within the example disclosed framework. Also discussed is the exact synthesis part using the framework introduced in Vadym Kliuchnikov and Jon Yard, “A framework for exact synthesis,” arXiv:1504.04350 (2015); U.S. Provisional Patent Application No. 62/146,182 filed on Apr. 10, 2015, and entitled “METHOD AND SYSTEM FOR QUANTUM CIRCUIT SYNTHESIS USING QUATERNION ALGEBRA”; and PCT International Application No. PCT/US2016/025958 filed on Apr. 5, 2016, and entitled “METHOD AND SYSTEM FOR QUANTUM CIRCUIT SYNTHESIS USING QUATERNION ALGEBRA”. This section follows the “Approximation methods based on exact synthesis” Section above (Section 2.2) and shows that Clifford+T can be described by totally definite quaternion algebra. Recall the following definition: Definition 2.3. a gate set

is described by quaternion algebra if the following data is defined:

TABLE 1 Summary of the results of running the algorithm for different gate sets. The quaternion algebra ${{{corresponding}\mspace{14mu}{to}\mspace{14mu}{the}\mspace{14mu}{gate}\mspace{14mu}{set}\mspace{14mu}{is}\mspace{14mu} Q} = {{\left( \frac{a,b}{F} \right)\mspace{14mu}{where}\mspace{14mu} F} = {{{\mathbb{Q}}\left( \zeta_{n} \right)}\bigcap{\mathbb{R}}}}};{N_{G}{\;\mspace{11mu}}{is}\mspace{14mu}{the}\mspace{14mu}{number}\mspace{14mu}{of}}$ generators to be used for exact synthesis; [L:L ∩ M] is the index of the intersection of L (the generalized Lipschitz order) and

 (the maximal order used to define the gate set) in L; N_(cnj) is the number of conjugacy classes of maximal orders in quaternion algebra Q; t_(exact) is the time in seconds spent on precomputation required for the exact synthesis part of the example algorithm; t_(approx) is the time in seconds spent on precomputation required for the approximation part; C_(min) is the additive constant appearing in Theorem 4.1; C_(min)/log N(

₁) is the ratio between C_(min) and the log of the norm of the ideal from S with the smallest norm (see Definition 2.3) GATE SET n a b N_(G) [L:L ∩

] N_(cnj) t_(approx) t_(exact) C_(min) C_(min)/log N(

₁) Clifford + T 8 −1 −1 3 1 1 0.160 1.360 6.099 8.799 V-basis 4 −1 −1 6 1 1 0.070 0.700 2.079 1.292 Clifford + R_(z) (π/8) 16 −1 −1 9 2 2 0.410 15.110 21.360 30.815 Clifford + T + V 8 −1 −1 29 1 1 0.170 5.680 6.099 8.799

(1) A quaternion gate set specification

F, σ, a, b,

, S

,

(2) A set of canonical generators

_(,S) of

_(S),

(3) A map Circuit as described by equation (1).

For Clifford+T gate set, one can choose:

={R _(α)(π/4),R _(α)(π/2):α=x,y,z}

It can separately be written R_(α)(π/2) because they generate Clifford group. Clifford gates are much cheaper in practice and typical cost function for Clifford+T gate set used in practice is: cos t(R _(α)(π/4))=1, cos t(R _(α)(π/2))=0, α=x,y,z.

Now consider the following quaternion gate set specification for Clifford+T:

-   -   F=         (ζ₈+ζ₈ ⁻¹) where ζ₈=e^(2iπ/8), let also θ be a primitive element         of F (in other words every element of F can be represented as         a₀+a₁θ where a₀,a₁ are rational numbers),     -   embedding σ:F→         is defined as σ(θ)=√{square root over (2)},     -   a=−1 and b=−1,     -   maximal order         of quaternion algebra

$Q = \left( \frac{{- 1},{- 1}}{F} \right)$

-   -    is

${{\mathbb{Z}}_{F} + {\frac{\left( {i + 1} \right)\theta}{2}{\mathbb{Z}}_{F}} + {\frac{\left( {j + 1} \right)\theta}{2}{\mathbb{Z}}_{F}} + {\frac{1 + i + j + k}{2}{\mathbb{Z}}_{F}}},$

-   -   where         _(F)=         [√{square root over (2)}] is a ring of integers of F,     -   S={         } where         =(2−θ)         _(F) (note that 2−θ is totally positive element of F).         Note that the discriminant of         is equal to         _(F), therefore         is coprime to it. This implies that the set         _(S) is infinite.

Using notation q_(z)=i, q_(y)=j, q_(x)=k, one can obtain set

_(Q) based on the following correspondence q _(t,α)=1+θ(1−q _(α))/2U _(q)(q _(t,α))=R _(α)(π/4) q _(c,α)=θ(1−q _(α))/2U _(q)(q _(c,α))=R _(α)(π/2) where α=x,y,z.

The next step is to compute

_(,S) using the algorithm from Vadym Kliuchnikov and Jon Yard, “A framework for exact synthesis,” arXiv:1504.04350 (2015). One can find that quaternion algebra Q has a trivial two sided ideal class group and that the number of conjugacy classes of maximal orders of Q is one. In this case, the situation is relatively straightforward. The set

_(,S) is equal to gen_(S)(

)∪gen_(u)(

). The set gen_(S)(

) comprises N(

)+1=3 elements with reduced norm a₀+a₁θ. The set gen_(u)(

) comprises three generators of the finite group of units of maximal order

modulo units of

_(F). The results of exemplary computations are the following:

-   -   gen_(S)(         )={q₁,q₂,q₃} where         q ₁=½+i/2+(θ−1)j/2+(θ−1)k/2,q ₂=(−θ+2)/2−θj/2,q         ₃=(−θ+2)j/2+θk/2.     -   gen_(u)(         )={u₁,u₂,u₃} where         u ₁=½−i/2−j/2−k/2,u ₂ =−j,u ₃=θ/2−θj/2.         Algorithmically, one can find that u₁=q_(c,z)q_(c,x), u₂=q_(c,y)         ², and u₃=q_(c,y). One sees that the unit group of         modulo units of         _(F) corresponds to Clifford group. Next, one finds that         q₁=q_(t,x)q_(c,z) ²q_(c,x)q_(c,z), q₂=−q_(t,z)q_(c,z)         ²q_(c,x)q_(c,z) and q₃=q_(t,y)q_(c,y) ³ up to a unit of         _(F). In general, the elements of set gen_(S)(         ) are defined up to right-hand side multiplication by a unit of         . For this reason, one can choose gen_(S)(         )={q_(t,α):α=x,y,z}. The map Circuit becomes almost trivial in         this case. In the next examples, this detail is omitted and         gen_(S)(         ) is written using generators convenient for the example         application.

Note that nrd q_(t,α)

_(F)=

therefore the cost vector corresponding to each q_(t,α) is (1). For all elements of the unit group the cost vector is (0). In this example original cost definition completely matches cost obtained based on cost vectors. Table 2, 5 shows the results of running the example circuit synthesis algorithm for Clifford+T gate set.

6.2. V-Basis.

In this section, a description of how to obtain results from Vadym Kliuchnikov and Jon Yard, “A framework for exact synthesis,” arXiv:1504.04350 (2015) within the example framework is provided.

TABLE 2 Results of running the algorithm for Clifford + T. Approximation of rotation R_(z)(0.1) with precision ε and cost vector (L₁). All columns except N_(tr,min) and N_(tr,max) are averages over 1000 runs of the algorithm; L′₁ is the averaged T-count of the found circuits; ρ(U_(q),R_(z)(ϕ)) is the obtained quality of approximation; N_(tr,min), N_(tr,max), N_(tr,avg) are minimum, maximum and average of the number of the main loop iterations in the procedure APPROXIMATE over all samples; t_(exact) and t_(approx) is time in seconds spent on approximation and exact synthesis stages of the algorithm. ε L₁ L′₁ ρ(U_(q), R_(z)(φ)) N_(tr,min) N_(tr,max) N_(tr,avg) t_(exact) t_(approx) 1 · 10⁻⁵  76  75.06 8.62 · 10⁻⁶ 1 101 19.21 0.023 0.073 1 · 10⁻⁶  89  87.88 8.65 · 10⁻⁷ 1  62 16.66 0.028 0.197 1 · 10⁻⁷ 102 100.88 8.66 · 10⁻⁸ 1 118 23.77 0.031 0.244 1 · 10⁻⁸ 116 114.94 8.67 · 10⁻⁹ 1 133 21.74 0.036 0.261 1 · 10⁻⁹ 129 127.94 8.66 · 10⁻¹⁰ 1 150 27.42 0.041 0.339 1 · 10⁻¹⁰ 142 141.06 8.64 · 10⁻¹¹ 1 264 31.32 0.047 0.386 1 · 10⁻¹⁵ 209 208.02  8.7 · 10⁻¹⁶ 2 193 49.1  0.077 0.507 1 · 10⁻²⁵ 341 339.88 8.67 · 10⁻²⁶ 1 508 87.9  0.154 1.219 1 · 10⁻³⁰ 408 407.04 8.67 · 10⁻³¹ 1 577 76.31 0.204 1.749

V-basis is defined using the following set

$G = \left\{ {\frac{1 \pm {2\;{iP}}}{\sqrt{5}},{{{iP}\text{:}\mspace{14mu} P} \in \left\{ {X,Y,Z} \right\}}} \right\}$

A typical cost function is:

${{{cost}\left( \frac{1 \pm {2\;{iP}}}{\sqrt{5}} \right)} = 1},{{{cost}({iP})} = 0},{{{for}\mspace{14mu} P} \in \left\{ {X,Y,Z} \right\}}$

The quaternion gate set specification is:

-   -   F=         ,     -   embedding σ:F→         is the only embedding of         into         ,     -   a=−1 and b=−1,     -   maximal order         of quaternion algebra

${Q = {{\left( \frac{{- 1},{- 1}}{\mathbb{Q}} \right)\mspace{14mu}{is}\mspace{14mu}{\mathbb{Z}}} + {i\;{\mathbb{Z}}} + {j\;{\mathbb{Z}}} + {\frac{1 + i + j + k}{2}{\mathbb{Z}}}}},$

-   -   S={         } where         =5         .         Note that the discriminant of         is equal to 2         and ideal         is co-prime to it. This implies that the set         _(S) is infinite.

Using notation q_(Z)=i, q_(Y)=j, q_(X)=k, one can obtain set

_(Q) based on the following correspondence: q _(V,±P)=1±2q _(P) U _(q)(q _(V,±P))=(1±2iP)/√{square root over (5)} U _(q)(q _(P))=iP where P∈{X,Y,Z}.

Similarly to the previous section, one can compute

_(,S) using the algorithm from Vadym Kliuchnikov and Jon Yard, “A framework for exact synthesis,” arXiv:1504.04350

TABLE 3 Results of running the algorithm for V basis. Approximation of rotation R_(z)(0.1) with precision ε and cost vector (L₁). All columns except N_(tr,min) and N_(tr,max) are averages over 1000 runs of the algorithm; L′₁ is the average V-count of the found circuits; ρ(U_(q), R_(z)(ϕ)) is the obtained quality of approximation; N_(tr,min), N_(tr,max), N_(tr,avg) are minimum, maximum and average of the number of the main loop iterations in the procedure APPROXIMATE over all samples; t_(exact) and t_(approx) is time in seconds spent on approximation and exact synthesis stages of the algorithm. ε L₁ L′₁ ρ(U_(q), R_(z)(φ)) N_(tr,min) N_(tr,max) N_(tr,avg) t_(exact) t_(approx) 1 · 10⁻⁵  30  30 8.62 · 10⁻⁶ 1 609 75.16 0.013 0.059 1 · 10⁻⁶  36  36 8.65 · 10⁻⁷ 1 622 94.88 0.017 0.074 1 · 10⁻⁷  42  42 8.66 · 10⁻⁸ 1 647 105.97 0.021 0.095 1 · 10⁻⁸  48  48 8.65 · 10⁻⁹ 1 808 123.1 0.024 0.109 1 · 10⁻⁹  53  53 8.64 · 10⁻¹⁰ 1 919 136.14 0.028 0.123 1 · 10⁻¹⁰  59  59 8.65 · 10⁻¹¹ 1 1,065 140.34 0.033 0.131 1 · 10⁻¹⁵  88  88 8.67 · 10⁻¹⁶ 1 1,470 227.13 0.058 0.225 1 · 10⁻²⁵ 145 145 8.65 · 10⁻²⁶ 1 3,064 380.11 0.129 0.435 1 · 10⁻³⁰ 173 173 8.62 · 10⁻³¹ 1 3,389 436.57 0.178 0.546 (2015). One can again find that quaternion algebra Q has a trivial two sided ideal class group and that the number of conjugacy classes of maximal orders of Q is one. The set

_(,S) is equal to gen_(S)(

)∪gen_(u)(

). The set gen_(S)(

) consists of N(

)+1=6 elements with reduced norm 5. The set gen_(u)(

) consists of two generators of the finite group of units of maximal order

modulo units of

. The results of the computations are the following:

-   -   gen_(S)(         )={q_(V,±P):P∈{X,Y,Z}}     -   gen_(u)(         )={i, (i+j+k+1)/2}

One of the generators of the unit group cannot be expressed as a product of elements of

_(Q). Indeed, all elements of

_(Q) belong to the Lipschitz order L=

+

i+

j+

k and (i+j+k+1)/2 does not. However, the example approximation algorithm finds q from L. It is possible to show that (in this particular example) the unit of

obtained in the end of exact synthesis of q must belong to L and therefore belongs to the subgroup of the unit group of

that contained in L. After a simple computation, one finds that this subgroup is generated by i,j,k.

Note that nrd q_(V,±P)

=

and the cost vector corresponding to each q_(V,±P) is (1). For all elements of the unit group the cost vector is (0). Similarly to Clifford+T case original cost definition completely matches cost obtained based on cost vectors. Table 3, 5 show the results of running the example circuit synthesis algorithm for V-basis.

6.3. Clifford+

$R_{z}\left( \frac{\pi}{8} \right)$ The approximation part of the result for this gate set is new. An exact synthesis algorithm for this gate set was first described in Simon Forest, David Gosset, Vadym Kliuchnikov, and David McKinnon, “Exact synthesis of single-qubit unitaries over clifford-cyclotomic gate sets,” Journal of Mathematical Physics, 56(8):082201 (2015) in the language of SO(3) representation of unitary matrices over the ring

[ζ₁₆, ½]. It can be shown that the output of the approximation stage of the example algorithm can be converted to a unitary matrix over

[ζ₁₆, ½]. Therefore the algorithm developed in Simon Forest, David Gosset, Vadym Kliuchnikov, and David McKinnon, “Exact synthesis of single-qubit unitaries over clifford-cyclotomic gate sets,” Journal of Mathematical Physics, 56(8):082201 (2015) can be applied instead of the exact synthesis algorithm for quaternions used here.

For Clifford+

$R_{z}\left( \frac{\pi}{8} \right)$ gate set, one can choose:

={R _(α)(±π/8),R _(α)(±3π/8),R _(α)(π/4),R _(α)(π/2):α=x,y,z}

For this example, the quaternion gate set specification is:

-   -   F=         (ζ₁₆+ζ₁₆ ⁻¹) where ζ₁₆=e^(2iπ/16), let also θ be a primitive         element of F (in other words every element of F can be         represented as a₀+ . . . +a₃θ where a₀, . . . , a₃ are rational         numbers),     -   embedding σ:F→         is defined as σ(θ)=2 cos(2π/16),     -   a=−1 and b=1,     -   maximal order         of quaternion algebra

$Q = \left( \frac{{- 1},{- 1}}{F} \right)$

${{\mathbb{Z}}_{F} + {\frac{\left( {i + 1} \right)\xi}{2}{\mathbb{Z}}_{F}} + {\frac{\left( {j + 1} \right)\xi}{2}{\mathbb{Z}}_{F}} + {\frac{1 + i + j + k}{2}{\mathbb{Z}}_{F}}},$

-   -   where ξ=θ²−2, σ(ξ)=√{square root over (2)} and         _(F)=         [2 cos(2π/16)] is a ring of integers of F,     -   S={         } where         =θ         _(F)=(−θ³+4θ²+θ−2)         _(F) (note that −θ³+4θ²+θ−2 is totally positive element of F).         The definition of maximal order         has essentially the same shape as the definition used for         Clifford+T case. The only difference is that it defined using         different ring of integers.

The discriminant of

is equal to

_(F), therefore

is coprime to it. This implies that the set

_(S) is infinite.

Using notation q_(z)=i, q_(y)=j, q_(x)=k, η=θ³−3θ (note σ(η)=2 sin(π/8)), one can obtain set

_(Q) based on the following correspondence: q _(t,α)=1+ξ(1−q _(α))/2U _(q)(q _(t,α))=R _(α)(π/4) q _(c,α)=(1−q _(α))/2U _(q)(q _(c,α))=R _(α)(π/2) q _(1/8,α)=θ(1+(θ−ηq _(α))/2)U _(q)(q _(1/8,α))=R _(α)(π/8) q _(3/8,α)=θ(1+(η−θq _(α))/2)U _(q)(q _(3/8,α))=R _(α)(3π/8) where α∈{x,y,z}.

The next step is to compute

_(,S) using the algorithm from Vadym Kliuchnikov and Jon Yard, “A framework for exact synthesis,” arXiv:1504.04350 (2015). One finds that quaternion algebra Q has a trivial two sided ideal class group and two different conjugacy classes of maximal orders of Q. The set

_(,S) is equal to gen_(S)(

)∪gen_(u)(

). The set gen_(u)(

) consists of three generators of the finite group of units of maximal order

modulo units of

. As one has two conjugacy classes of maximal orders, it desirable to build an ideal principality graph (which is a tree in this case, see FIG. 16) to find the set gen_(S)(

). The result of the computation is the following:

-   -   gen_(S)(         )={q_(t,α),q_(1/8,α),q_(3/8,α):α∈{x,y,z} }     -   gen_(u)(         )={q_(c,α):α∈{x,y,z}}

In particular, FIG. 16 is an ideal principality graph 1600 used to find generators for exact synthesis algorithm. The illustrated tree corresponds to quaternion gate set specification for Clifford+R_(z)(π/8) gate set. See also Vadym Kliuchnikov and Jon Yard, “A framework for exact synthesis,” arXiv:1504.04350 (2015)

This computation reproduces the result from Simon Forest, David Gosset, Vadym Kliuchnikov, and David McKinnon, “Exact synthesis of single-qubit unitaries over clifford-cyclotomic gate sets.” Journal of Mathematical Physics, 56(8):082201 (2015) showing that all matrices over the ring

[ζ₁₆, ½] can be exactly represented using gate set

. Because there are two conjugacy classes of maximal orders, the situation with the cost of generators becomes more interesting. For quaternions q_(t,α), nrd(q_(t,α))=

² and their cost vector is (2). For other elements of gen_(S)(

), nrd(q_(1/8,α))=

³ and nrd(q_(3/8,α))=

³ and their cost vector is (3). In the case when S contains only one prime ideal, the cost of each generator from gen_(u)(

) is precisely equal to the distance from the root to corresponding node. Above cost values also reproduce results from Simon Forest, David Gosset, Vadym Kliuchnikov, and David McKinnon, “Exact synthesis of single-qubit unitaries over clifford-cyclotomic gate sets,” Journal of Mathematical Physics, 56(8):082201 (2015).

Note that while approximating, one only has control over the overall value of the cost vector. If one requested cost L, then the result can have any number L_(t) of the T gates and any number L_(1/8,3/8) of R_(z)(π/8) and R_(z)(3π/8) rotations as soon as L=2L_(t)+3L_(1/8,3/8). As usual, the cost of Clifford gates is assumed to be zero.

Another interesting aspect of this example is that the generalized Lipschitz order is not contained in maximal order M above. Orders in a totally definite quaternion algebra can be given a structure of the lattice using bilinear form

(q₁q₂*). One can find that the index of sub-lattice L∩M in L is two. This means that half of the points from L belongs L∩M. In the example approximation algorithm, the result is tested to determine if it is in M in the end. If this is not the case, the procedure tries again. Experiments show that the procedure produces a quaternion in L∩M in half of the experiments. Tables 4, 5 show the results of running the example circuit synthesis algorithm for Clifford+R_(z)(π/8) gate set.

6.4. Clifford+T+V.

For Clifford+T+V gate set, one can choose set G to be:

={R _(α)(π/4),R _(α)(±2a tan(2)),R _(α)(π/2):α=x,y,z} It is not difficult to check that R_(α)(±2a tan(2)) correspond to 6 V gates.

TABLE 4 Results of running the example algorithm for Clifford + R_(z)(π/8). Part 1. Approximation of rotation R_(z)(0.1) with precision ε and cost vector (L₁). All columns except N_(tr,min) and N_(tr,max) are averages over 1000 runs of the algorithm; (L′₁) is a cost vector of a found circuit; ρ(U_(q), R_(z)(ϕ)) is the obtained quality of approximation; N_(tr,min), N_(tr,max), N_(tr,avg) are minimum, maximum and average of the number of the main loop iterations in the procedure APPROXIMATE over all samples; t_(exact) and t_(approx) is time in seconds spent on approximation and exact synthesis stages of the algorithm. ε L₁ L′₁ ρ(U_(q), R_(z)(φ)) N_(tr,min) N_(tr,max) N_(tr,avg) t_(exact) t_(approx) 1 · 10⁻⁵  98  96 8.65 · 10⁻⁶ 1 151 19.56 0.064 0.860 1 · 10⁻⁶ 111 108.98 8.66 · 10⁻⁷ 1 171 21.94 0.075 0.951 1 · 10⁻⁷ 124 121.95 8.66 · 10⁻⁸ 1 168 22.69 0.085 1.049 1 · 10⁻⁸ 138 136.03 8.67 · 10⁻⁹ 1 217 24.59 0.098 1.207 1 · 10⁻⁹ 151 148.96 8.66 · 10⁻¹⁰ 1 207 27.75 0.108 1.349 1 · 10⁻¹⁰ 164 161.92 8.65 · 10⁻¹¹ 1 223 29.62 0.118 1.541 1 · 10⁻¹⁵ 231 229.01 8.66 · 10⁻¹⁶ 1 267 40.74 0.168 2.057 1 · 10⁻²⁵ 364 361.93 8.66 · 10⁻²⁶ 1 422 59.35 0.273 4.647 1 · 10⁻³⁰ 430 428.04 8.66 · 10⁻³¹ 1 466 68.1 0.329 6.551

TABLE 5 Results of running the example algorithm for Clifford + R_(z)(π/8). Part 2. Approximation of rotation R_(z)(0.1) with precision ε and cost vector (L₁). All columns except N_(L∩M), L′_(t), L′_(1/8,3/8) are averages over 1000 runs of the algorithm; N_(L∩M) is the number of outputs of the procedure APPROXIMATE that are in the maximal order  

 ; (L′₁) is a cost vector of a found circuit; ρ(U_(q), R_(z)(ϕ)) is the obtained quality of approximation; L′_(t) is the number of T gates in the resulting circuit (averaged over outputs of the procedure APPROXIMATE that are in the maximal order  

 ); L′_(1/8,3/8) is the number of R_(z)(π/8) and R_(z)(3π/8) gates in the resulting circuit (averaged in the same way as L_(t)′). ε L₁ L′₁ ρ(U_(q), R_(z)(φ)) L′_(t) L′_(1/8,3/8) N_(L∩M) 1 · 10⁻⁵  98  96 8.65 · 10⁻⁶ 19.49 19.01 509 1 · 10⁻⁶ 111 108.98 8.66 · 10⁻⁷ 22.48 21.34 518 1 · 10⁻⁷ 124 121.95 8.66 · 10⁻⁸ 25.26 23.81 478 1 · 10⁻⁸ 138 136.03 8.67 · 10⁻⁹ 28.09 26.62 484 1 · 10⁻⁹ 151 148.96 8.66 · 10⁻¹⁰ 30.83 29.1 480 1 · 10⁻¹⁰ 164 161.92 8.65 · 10⁻¹¹ 32.74 32.15 517 1 · 10⁻¹⁵ 231 229.01 8.66 · 10⁻¹⁶ 46.38 45.42 491 1 · 10⁻²⁵ 364 361.93 8.66 · 10⁻²⁶ 73.27 71.8 474 1 · 10⁻³⁰ 430 428.04 8.66 · 10⁻³¹ 85.94 85.39 528

The gate set specification is similar to Clifford+T case except for the set S.

-   -   F=         (ζ₈+ζ₈ ⁻¹) where ζ₈=e^(2iπ/8), let also θ be a primitive element         of F (in other words every element of F can be represented as         a₀+a₁θ where a₀,a₁ are rational numbers),     -   embedding σ:F→         is defined as σ(θ)=√{square root over (2)},     -   a=−1 and b=−1,     -   maximal order         of quaternion algebra

$Q = \left( \frac{{- 1},{- 1}}{F} \right)$ is

${{\mathbb{Z}}_{F} + {\frac{\left( {i + 1} \right)\theta}{2}{\mathbb{Z}}_{F}} + {\frac{\left( {j + 1} \right)\theta}{2}{\mathbb{Z}}_{F}} + {\frac{1 + i + j + k}{2}{\mathbb{Z}}_{F}}},$

-   -   where         =         [√{square root over (2)}] is a ring of integers of F,     -   S={         ₁,         ₂} where         ₁=(2−θ)         _(F) and         ₂=5         _(F).

Using notation q_(z)=i, q_(y)=j, q_(x)=k, one can obtain set

_(Q) based on the following correspondence q _(t,α)=1+θ(1−q _(α))/2U _(q)(q _(t,α))=R _(α)(π/4) q _(v,±α)=1∓2q _(α) U _(q)(q _(v,±α))=R _(α)(±2a tan(2)) q _(c,α)=θ(1−q _(α))/2U _(q)(q _(c,α))=R _(α)(π/2) where α=x, y, z.

The next step is to compute

_(,S) using the algorithm from Vadym Kliuchnikov and Jon Yard, “A framework for exact synthesis,” arXiv:1504.04350 (2015). One can find that quaternion algebra Q has a trivial two sided ideal class group and that the number of conjugacy classes of maximal orders of Q is one. The set

_(,S) is equal to gen_(S)(

)∪gen_(u)(

). The set gen_(S)(

) consists of N(

₁)+1=3 elements with reduced norm 2−θ and N(

₂)+1=26 elements with reduced norm 5. The set gen_(u)(

) consists of three generators of the finite group of units of maximal order

modulo units of

_(F) and is the same as in Clifford+T case because maximal order

is the same.

Here, the set gen_(S)(

) is discussed in more detail. One can find that it contains quaternions corresponding to R_(α)(π/4) gates and all 6 V gates. One is then left with 20 quaternions with reduced norm 5 that were not in the set

. These can be expressed in terms of elements of

. The following equivalence relation on quaternions can then be introduced: q ₁ ˜q ₂ if and only if q ₁ =u ₁ q ₂ u ₂ for u ₁ ,u ₂−units of M In this case, it means that corresponding unitaries are equivalent up to a Clifford and therefore will have the same cost of implementation. There are four equivalence classes in gen_(S)(

) corresponding to the relation ˜. Two of them are {q_(t,α):α∈{x,y,z}} and {q_(v,±α):α∈{x,y,z}}. The remaining twenty quaternions with reduced norm 5 split into two classes c₁ and c₂ of size 8 and 12. Next, one can find that all quaternions from c₂ are equal to u ₁ q _(t,α(1)) q _(v,±α(2)) u ₂ q _(t,α(3)) ⁻¹ where u ₁ ,u ₂ are units of

,α(k)∈{x,y,z},k=1,2,3.

The quaternions from the set c₁ can be expressed as u ₁ q _(t,α(1)) q _(t,α(2)) q _(v,±α(3)) u ₂ q _(t,α(4)) ⁻¹ q _(t,α(5)) ⁻¹ where u ₁ ,u ₂ are units of

,α(k)∈{x,y,z}k=1, . . . ,5.

TABLE 6 Results of running the example algorithm for Clifford + T + V. Approximation of rotation R_(z)(0.1) with precision ε and cost vector (L₁, L₂). All columns except N_(tr,min) and N_(tr,max) are averages over 1000 runs of the algorithm; (L′₁, L′₂) is a cost vector of a found circuit; ρ(U_(q), R_(z)(ϕ)) is the obtained quality of approximation; N_(tr,min), N_(tr,max), N_(tr,avg) are minimum, maximum and average of the number of the main loop iterations in the procedure APPROXIMATE over all samples; t_(exact) and t_(approx) is time in seconds spent on approximation and exact synthesis stages of the algorithm. ε L₁ L₂ L′₁ L′₂ ρ(U_(q),R_(z)(φ)) N_(tr,min) N_(tr,max) N_(tr,avg) t_(exact) t_(approx) 1 · 10⁻⁵  30 10  29.06 10 8.65 · 10⁻⁶ 1 149 18.31 0.075 0.070 1 · 10⁻⁶  36 12  35.02 12 8.64 · 10⁻⁷ 1 158 25.07 0.089 0.241 1 · 10⁻⁷  42 14  40.98 14 8.65 · 10⁻⁸ 1 204 28.59 0.107 0.277 1 · 10⁻⁸  48 16  47 16 8.66 · 10⁻⁹ 1 210 32.97 0.129 0.366 1 · 10⁻⁹  51 17  50.02 17 8.65 · 10⁻¹⁰ 1 223 29.32 0.140 0.367 1 · 10⁻¹⁰  57 19  55.99 19 8.66 · 10⁻¹¹ 1 347 37.59 0.155 0.429 1 · 10⁻¹⁵  84 28  83.03 28 8.66 · 10⁻¹⁶ 1 383 52.66 0.233 0.559 1 · 10⁻²⁵ 135 45 133.96 45 8.65 · 10⁻²⁶ 1 702 83.11 0.408 1.306 1 · 10⁻³⁰ 162 54 161.05 54 8.66 · 10⁻³¹ 1 618 101.72 0.508 1.978 In practice it can be more beneficial to design circuits for all 26 gates corresponding to quaternions with norm 5 directly, because T gates are usually expensive to implement. Table 6 shows the results of running the example circuit synthesis algorithm for Clifford+T+V gate set.

7. Other Examples of Using Approximation Algorithm

In this section, results of running an example implementation of the disclosed approximation algorithm for a series of quaternion gate sets specification are provided. Recall that the quaternion gate set specification is:

Definition 2.2.

A quaternion gate set specification is a tuple

F, σ, a, b,

, S

where:

-   -   F is a totally real number field and σ is an embedding of F into     -   a, b are elements of F that define the quaternion algebra

$\left( \frac{a,b}{F} \right)$

-   -    over F     -   is a maximal order of

$\left( \frac{a,b}{F} \right)$

-   -   S={         ₁, . . . ,         _(M)} is a set of prime ideals of F

The family of examples is parametrized by even integer n. Number field F corresponds to the real subfield of cyclotomic field

(ζ_(n)+ζ_(n) ⁻¹) with primitive element θ. Element a of F is chosen such that relative extension F(a) is a cyclotomic field, b=−1. The approximation part of the example algorithm is independent on maximal order

, so there was no restriction to any specific choice of

. Set S contains one prime ideal above 2. If there is more than one such ideal, it was chosen at random. The tables below summarize the series of examples using an implementation of the disclosed algorithm.

TABLE 7 n a b t_(approx) C_(min) C_(min)/log N (p₁) 8 −1 −1 0.150 6.099 8.799 10  θ-3 −1 0.080 4.812 3.471 12 −1 −1 0.100 7.100 10.244 14 θ²-4 −1 0.140 9.590 4.612 16 −1 −1 0.430 21.360 30.815 18 θ²-4 −1 0.130 10.679 5.136 20 −1 −1 0.260 20.598 14.859 22 θ²-4 −1 0.680 24.440 7.052 24 −1 −1 0.290 21.766 31.402 26 θ²-4 −1 0.460 34.891 8.389 28 −1 −1 0.540 40.915 19.676 30 θ²-4 −1 0.320 17.027 6.141 32 −1 −1 1.180 71.573 103.259 34 θ²-4 −1 1.060 59.243 21.367 36 −1 −1 0.440 41.956 20.177 38 θ²-4 −1 1.350 72.404 11.606 40 −1 −1 1.250 67.612 48.772 42 θ²-4 −1 0.480 35.413 8.515 44 −1 −1 2.730 97.529 28.141

This table (and the tables below) summarize the results of running offline part of the example algorithm: n is the number of example in the family described above; t_(approx) is the time in seconds spent on the offline stage required for the approximation part; C_(min) is the additive constant appearing in Theorem 4.1; C/log N(

₁) is the ratio between C_(min) and the log of the norm of the ideal in S. Next, tables are shown with the averages over 100 runs of the example algorithm with different target precisions ε and target cost vector (L₁) and target angle φ=0.1 for each example for n=8, 10, . . . , 44. All columns of the tables except N_(tr,min) and N_(tr,max) are averages over 100 runs of the algorithm; ρ(U_(q),R_(z)(φ)) is the obtained quality of approximation; N_(tr,min),N_(tr,max),N_(tr,avg) are minimum, maximum and average of the number of the main loop iterations in the procedure APPROXIMATE over all samples; t_(approx) is time in seconds spent on online part of the approximation stage of the algorithm.

TABLE 8 ε L₁ ρ(U_(q), R_(z)((φ)) N_(tr,min) N_(tr,max) N_(tr,avg) t_(approx) n = 8 1 · 10⁻⁵  76 8.65 · 10⁻⁶ 1 128 19.03 0.073 1 · 10⁻⁶  89 8.69 · 10⁻⁷ 1 112 26.29 0.212 1 · 10⁻⁷ 102 8.66 · 10⁻⁸ 1 163 29.55 0.242 1 · 10⁻⁸ 116 8.63 · 10⁻⁹ 1 151 31.67 0.271 1 · 10⁻⁹ 129 8.62 · 10⁻¹⁰ 2 195 31.49 0.327 1 · 10⁻¹⁰ 142 8.67 · 10⁻¹¹ 1 186 36.84 0.376 1 · 10⁻¹⁵ 209 8.66 · 10⁻¹⁶ 1 299 52.42 0.478 1 · 10⁻²⁵ 341 8.68 · 10⁻²⁶ 2 354 78.44 1.111 1 · 10⁻³⁰ 408 8.69 · 10⁻³¹ 1 397 98.67 1.750 n = 10 1 · 10⁻⁵  37 8.61 · 10⁻⁶ 1 102 19.43 0.099 1 · 10⁻⁶  44 8.62 · 10⁻⁷ 1 160 42.21 0.277 1 · 10⁻⁷  50 8.68 · 10⁻⁸ 1 271 49.29 0.389 1 · 10⁻⁸  57 8.64 · 10⁻⁹ 1 139 25.76 0.350 1 · 10⁻⁹  64 8.68 · 10⁻¹⁰ 2 292 62.61 0.558 1 · 10⁻¹⁰  70 8.62 · 10⁻¹¹ 1 312 71.45 0.681 1 · 10⁻¹⁵ 104 8.64 · 10⁻¹⁶ 2 570 108.9 0.935 1 · 10⁻²⁵ 170 8.67 · 10⁻²⁶ 1 1,321 185.05 2.295 1 · 10⁻³⁰ 203 8.64 · 10⁻³¹ 1 509 90.78 2.589 n = 12 1 · 10⁻⁵  78 8.68 · 10⁻⁶ 1 67 17.96 0.069 1 · 10⁻⁶  90 8.65 · 10⁻⁷ 1 76 19.69 0.200 1 · 10⁻⁷ 104 8.66 · 10⁻⁸ 1 105 27.02 0.233 1 · 10⁻⁸ 118 8.68 · 10⁻⁹ 1 174 26.97 0.291 1 · 10⁻⁹ 130 8.61 · 10⁻¹⁰ 1 146 27.91 0.317 1 · 10⁻¹⁰ 144 8.67 · 10⁻¹¹ 1 149 33.79 0.371 1 · 10⁻¹⁵ 210 8.59 · 10⁻¹⁶ 1 238 46.1 0.461 1 · 10⁻²⁵ 344 8.66 · 10⁻²⁶ 1 351 71.82 1.164 1 · 10⁻³⁰ 410 8.64 · 10⁻³¹ 2 527 87.47 1.633 n = 14 1 · 10⁻⁵  27 8.65 · 10⁻⁶ 1 87 20.43 0.163 1 · 10⁻⁶  32 8.63 · 10⁻⁷ 1 206 29.84 0.497 1 · 10⁻⁷  36 8.67 · 10⁻⁸ 1 157 31.2 0.538 1 · 10⁻⁸  41 8.67 · 10⁻⁹ 1 198 36.99 0.708 1 · 10⁻⁹  45 8.65 · 10⁻¹⁰ 1 256 37.85 0.783 1 · 10⁻¹⁰  49 8.62 · 10⁻¹¹ 1 199 39.5 0.807 1 · 10⁻¹⁵  72 8.66 · 10⁻¹⁶ 1 235 51.54 1.107 1 · 10⁻²⁵ 116 8.66 · 10⁻²⁶ 1 638 85.51 2.615 1 · 10⁻³⁰ 138 8.64 · 10⁻³¹ 5 459 111.86 4.064

TABLE 9 ε L₁ ρ(U_(q), R_(z)((φ)) N_(tr,min) N_(tr,max) N_(tr,avg) t_(approx) n = 16 1 · 10⁻⁵  98 8.66 · 10⁻⁶ 1 110 22.16 0.816 1 · 10⁻⁶ 111 8.66 · 10⁻⁷ 1 99 19.03 0.911 1 · 10⁻⁷ 124 8.67 · 10⁻⁸ 1 116 22.55 0.987 1 · 10⁻⁸ 138 8.66 · 10⁻⁹ 1 126 22.24 1.077 1 · 10⁻⁹ 151 8.67 · 10⁻¹⁰ 1 86 24.06 1.192 1 · 10⁻¹⁰ 164 8.66 · 10⁻¹¹ 1 241 38.3 1.504 1 · 10⁻¹⁵ 231 8.67 · 10⁻¹⁶ 1 221 34.33 1.856 1 · 10⁻²⁵ 364 8.69 · 10⁻²⁶ 2 352 58.02 4.336 1 · 10⁻³⁰ 430 8.67 · 10⁻³¹ 1 456 69.24 6.612 n = 18 1 · 10⁻⁵  28 8.67 · 10⁻⁶ 1 77 22.13 0.379 1 · 10⁻⁶  32 8.63 · 10⁻⁷ 1 133 34.03 0.445 1 · 10⁻⁷  37 8.61 · 10⁻⁸ 1 45 12.65 0.549 1 · 10⁻⁸  41 8.69 · 10⁻⁹ 1 81 16.9 0.667 1 · 10⁻⁹  45 8.65 · 10⁻¹⁰ 1 103 16.36 0.653 1 · 10⁻¹⁰  50 8.67 · 10⁻¹¹ 3 162 39.05 0.720 1 · 10⁻¹⁵  72 8.67 · 10⁻¹⁶ 1 358 49.58 0.944 1 · 10⁻²⁵ 116 8.66 · 10⁻²⁶ 1 476 82.59 2.280 1 · 10⁻³⁰ 139 8.66 · 10⁻³¹ 1 226 45.13 3.827 n = 20 1 · 10⁻⁵  50 8.64 · 10⁻⁶ 1 49 12.23 0.833 1 · 10⁻⁶  56 8.62 · 10⁻⁷ 1 37 10.13 0.876 1 · 10⁻⁷  62 8.65 · 10⁻⁸ 1 72 14.3 0.974 1 · 10⁻⁸  70 8.63 · 10⁻⁹ 1 76 15.4 1.151 1 · 10⁻⁹  76 8.66 · 10⁻¹⁰ 1 82 16.41 1.225 1 · 10⁻¹⁰  82 8.68 · 10⁻¹¹ 1 176 18.24 1.366 1 · 10⁻¹⁵ 116 8.66 · 10⁻¹⁶ 1 124 21.75 1.916 1 · 10⁻²⁵ 182 8.66 · 10⁻²⁶ 1 154 26.96 4.207 1 · 10⁻³⁰ 216 8.65 · 10⁻³¹ 1 352 41.98 6.611 n = 22 1 · 10⁻⁵  21 8.67 · 10⁻⁶ 1 51 13.38 1.544 1 · 10⁻⁶  23 8.69 · 10⁻⁷ 1 71 13.04 1.634 1 · 10⁻⁷  26 8.65 · 10⁻⁸ 1 78 22.48 1.863 1 · 10⁻⁸  29 8.66 · 10⁻⁹ 1 148 21.94 2.321 1 · 10⁻⁹  31 8.66 · 10⁻¹⁰ 2 77 28.86 2.612 1 · 10⁻¹⁰  34 8.66 · 10⁻¹¹ 1 105 31.44 2.423 1 · 10⁻¹⁵  47 8.65 · 10⁻¹⁶ 1 117 30.2 3.692 1 · 10⁻²⁵  74 8.64 · 10⁻²⁶ 1 246 56.32 8.760 1 · 10⁻³⁰  87 8.64 · 10⁻³¹ 1 329 65.86 12.830

TABLE 10 ε L₁ ρ(U_(q), R_(z)((φ)) N_(tr,min) N_(tr,max) N_(tr,avg) t_(approx) n = 24 1 · 10⁻⁵  98 8.65 · 10⁻⁶ 1 50 11.88 0.815 1 · 10⁻⁶ 112 8.66 · 10⁻⁷ 1 82 14.89 0.922 1 · 10⁻⁷ 126 8.68 · 10⁻⁸ 1 61 13.42 1.046 1 · 10⁻⁸ 138 8.65 · 10⁻⁹ 1 91 17.72 1.155 1 · 10⁻⁹ 152 8.68 · 10⁻¹⁰ 1 120 19.03 1.340 1 · 10⁻¹⁰ 166 8.66 · 10⁻¹¹ 1 159 23.67 1.240 1 · 10⁻¹⁵ 232 8.65 · 10⁻¹⁶ 1 146 23.65 1.919 1 · 10⁻²⁵ 364 8.64 · 10⁻²⁶ 2 192 37.63 4.528 1 · 10⁻³⁰ 432 8.67 · 10⁻³¹ 1 287 47.19 6.778 n = 26 1 · 10⁻⁵  20 8.66 · 10⁻⁶ 1 116 26.32 2.984 1 · 10⁻⁶  22 8.66 · 10⁻⁷ 1 153 35.34 3.472 1 · 10⁻⁷  24 8.66 · 10⁻⁸ 1 133 36.1 3.787 1 · 10⁻⁸  27 8.65 · 10⁻⁹ 1 135 26.2 3.539 1 · 10⁻⁹  29 8.67 · 10⁻¹⁰ 1 160 29.87 4.035 1 · 10⁻¹⁰  31 8.66 · 10⁻¹¹ 2 130 36.41 4.600 1 · 10⁻¹⁵  42 8.67 · 10⁻¹⁶ 1 291 55.91 7.531 1 · 10⁻²⁵  64 8.66 · 10⁻²⁶ 1 454 73.41 15.297 1 · 10⁻³⁰  75 8.66 · 10⁻³¹ 1 454 56.26 19.222 n = 28 1 · 10⁻⁵  42 8.66 · 10⁻⁶ 1 76 17.85 2.784 1 · 10⁻⁶  47 8.66 · 10⁻⁷ 1 66 14.99 2.548 1 · 10⁻⁷  51 8.68 · 10⁻⁸ 1 85 22.34 2.846 1 · 10⁻⁸  56 8.66 · 10⁻⁹ 1 128 20.5 3.168 1 · 10⁻⁹  60 8.67 · 10⁻¹⁰ 1 162 23.79 3.299 1 · 10⁻¹⁰  64 8.66 · 10⁻¹¹ 1 135 31.42 3.772 1 · 10⁻¹⁵  87 8.66 · 10⁻¹⁶ 1 224 31.28 5.903 1 · 10⁻²⁵ 131 8.66 · 10⁻²⁶ 1 204 43.22 12.796 1 · 10⁻³⁰ 153 8.65 · 10⁻³¹ 1 263 45.78 17.147 n = 30 1 · 10⁻⁵  23 8.65 · 10⁻⁶ 1 67 11.67 0.625 1 · 10⁻⁶  27 8.67 · 10⁻⁷ 1 61 14.29 0.825 1 · 10⁻⁷  30  8.7 · 10⁻⁸ 1 75 13.18 0.874 1 · 10⁻⁸  33 8.63 · 10⁻⁹ 1 63 16.27 0.982 1 · 10⁻⁹  37 8.68 · 10⁻¹⁰ 1 70 19.81 1.189 1 · 10⁻¹⁰  40 8.64 · 10⁻¹¹ 1 156 20.87 1.250 1 · 10⁻¹⁵  56 8.67 · 10⁻¹⁶ 1 154 31.02 1.574 1 · 10⁻²⁵  90 8.66 · 10⁻²⁶ 1 197 36.21 3.900 1 · 10⁻³⁰ 106 8.65 · 10⁻³¹ 1 210 49.94 5.623

TABLE 11 ε L₁ ρ(U_(q), R_(z)((φ)) N_(tr,min) N_(tr,max) N_(tr,avg) t_(approx) n = 32 1 · 10⁻⁵ 170 8.66 · 10⁻⁶ 1 133 30.76 8.326 1 · 10⁻⁶ 183 8.66 · 10⁻⁷ 1 182 35.13 8.556 1 · 10⁻⁷ 197 8.66 · 10⁻⁸ 1 415 36.98 10.075 1 · 10⁻⁸ 210 8.66 · 10⁻⁹ 1 191 31.61 10.317 1 · 10⁻⁹ 223 8.66 · 10⁻¹⁰ 1 178 41.23 11.530 1 · 10⁻¹⁰ 237 8.66 · 10⁻¹¹ 4 237 48.6 12.149 1 · 10⁻¹⁵ 303 8.66 · 10⁻¹⁶ 1 198 40.82 16.668 1 · 10⁻²⁵ 436 8.66 · 10⁻²⁶ 1 302 61.2 32.528 1 · 10⁻³⁰ 502 8.66 · 10⁻³¹ 1 346 69.61 42.645 n = 34 1 · 10⁻⁵  38 8.65 · 10⁻⁶ 1 238 33.2 7.857 1 · 10⁻⁶  42 8.66 · 10⁻⁷ 1 251 39.74 9.701 1 · 10⁻⁷  45 8.66 · 10⁻⁸ 1 162 36.29 9.778 1 · 10⁻⁸  48 8.65 · 10⁻⁹ 1 133 34.7 9.225 1 · 10⁻⁹  52 8.66 · 10⁻¹⁰ 1 206 37.21 10.723 1 · 10⁻¹⁰  55 8.66 · 10⁻¹¹ 1 198 45.27 13.598 1 · 10⁻¹⁵  72 8.66 · 10⁻¹⁶ 1 441 46.66 18.113 1 · 10⁻²⁵ 105 8.66 · 10⁻²⁶ 1 254 54.74 35.632 1 · 10⁻³⁰ 122 8.66 · 10⁻³¹ 1 394 77.65 51.990 n = 36 1 · 10⁻⁵  44 8.65 · 10⁻⁶ 1 109 19.75 2.746 1 · 10⁻⁶  48 8.65 · 10⁻⁷ 1 72 19.58 2.746 1 · 10⁻⁷  52 8.66 · 10⁻⁸ 1 105 23.99 2.879 1 · 10⁻⁸  56 8.66 · 10⁻⁹ 1 60 19.23 3.048 1 · 10⁻⁹  62 8.66 · 10⁻¹⁰ 1 152 30.17 3.880 1 · 10⁻¹⁰  66 8.66 · 10⁻¹¹ 1 160 28.35 4.004 1 · 10⁻¹⁵  88 8.67 · 10⁻¹⁶ 1 211 38.6 6.212 1 · 10⁻²⁵ 132 8.66 · 10⁻²⁶ 1 241 48.72 12.719 1 · 10⁻³⁰ 154 8.66 · 10⁻³¹ 1 311 61.65 18.355 n = 38 1 · 10⁻⁵  19 8.66 · 10⁻⁶ 1 149 26.57 13.018 1 · 10⁻⁶  21 8.66 · 10⁻⁷ 1 136 27.4 15.025 1 · 10⁻⁷  22 8.66 · 10⁻⁸ 1 208 39.52 15.637 1 · 10⁻⁸  24 8.66 · 10⁻⁹ 3 191 38.08 17.399 1 · 10⁻⁹  25 8.66 · 10⁻¹⁰ 1 157 35.81 18.434 1 · 10⁻¹⁰  27 8.66 · 10⁻¹¹ 1 180 39.61 19.639 1 · 10⁻¹⁵  34 8.66 · 10⁻¹⁶ 1 327 54.22 27.450 1 · 10⁻²⁵  49 8.67 · 10⁻²⁶ 2 361 65.16 56.848 1 · 10⁻³⁰  56 8.66 · 10⁻³¹ 1 459 71.79 73.107

TABLE 12 ε L₁ ρ(U_(q), R_(z)((φ)) N_(tr,min) N_(tr,max) N_(tr,avg) t_(approx) n = 40 1 · 10⁻⁵  82 8.66 · 10⁻⁶ 1 82 16.13 8.625 1 · 10⁻⁶  90 8.65 · 10⁻⁷ 1 79 18.79 9.455 1 · 10⁻⁷  96 8.66 · 10⁻⁸ 1 89 16.4 9.899 1 · 10⁻⁸ 102 8.66 · 10⁻⁹ 1 79 19.49 10.597 1 · 10⁻⁹ 110 8.66 · 10⁻¹⁰ 1 107 19.77 12.147 1 · 10⁻¹⁰ 116 8.67 · 10⁻¹¹ 1 131 16.77 11.968 1 · 10⁻¹⁵ 150 8.66 · 10⁻¹⁶ 1 130 25.94 18.669 1 · 10⁻²⁵ 216 8.66 · 10⁻²⁶ 1 181 35.68 33.633 1 · 10⁻³⁰ 250 8.66 · 10⁻³¹ 1 153 36.82 48.701 n = 42 1 · 10⁻⁵  20 8.66 · 10⁻⁶ 1 82 11.88 2.396 1 · 10⁻⁶  22 8.66 · 10⁻⁷ 1 64 14.22 2.547 1 · 10⁻⁷  25 8.65 · 10⁻⁸ 1 73 17.03 2.646 1 · 10⁻⁸  27 8.66 · 10⁻⁹ 1 101 16.12 2.950 1 · 10⁻⁹  29 8.67 · 10⁻¹⁰ 1 68 15.15 3.109 1 · 10⁻¹⁰  31 8.66 · 10⁻¹¹ 1 75 18.68 3.420 1 · 10⁻¹⁵  42 8.68 · 10⁻¹⁶ 1 173 22.95 5.191 1 · 10⁻²⁵  64 8.66 · 10⁻²⁶ 1 131 33.37 11.135 1 · 10⁻³⁰  75 8.66 · 10⁻³¹ 1 223 35.94 15.043 n = 44 1 · 10⁻⁵  42 8.66 · 10⁻⁶ 1 150 22.14 22.552 1 · 10⁻⁶  46 8.66 · 10⁻⁷ 1 141 25.34 25.966 1 · 10⁻⁷  48 8.66 · 10⁻⁸ 1 104 24.45 27.323 1 · 10⁻⁸  50 8.66 · 10⁻⁹ 1 112 29.35 27.838 1 · 10⁻⁹  54 8.66 · 10⁻¹⁰ 1 141 24.78 31.356 1 · 10⁻¹⁰  56 8.66 · 10⁻¹¹ 1 122 28.01 33.021 1 · 10⁻¹⁵  70 8.66 · 10⁻¹⁶ 1 163 27.98 44.929 1 · 10⁻²⁵  96 8.66 · 10⁻²⁶ 1 207 41 81.426 1 · 10⁻³⁰ 108 8.66 · 10⁻³¹ 2 141 38.83 105.410

8. Further Implementation Details

In FIG. 11 and FIG. 12, pseudo-code is provided for the versions of procedures UNIT-ANDJUST (FIG. 11) and RANDOM-INTEGER-POINT (FIG. 12) that can be used in implementations of the disclosed algorithm. They are useful to reproduce values of additive constant listed in Section 6 and Section 7. They both rely on the Nearest Plane Algorithm shown in FIG. 10). For given target vector t, the Nearest Plane algorithm finds unique lattice vector Bv inside t+C(B*). See Section 3.5 for the definition of C(B*) and other related definitions.

9. Example Overall Quantum Computer Gate Set Generation and Implementation Techniques

FIG. 18 is a flow chart 1800 showing an example synthesis process for producing a gate set in accordance with embodiments of the disclosed technology. The process can be performed, for instance, by a quantum computer synthesis tool adapted for use in a quantum computer design and implementation process and implemented by one or more computing devices. Although the process in FIG. 18 is described with respect to a single qubit gate sets, the size and complexity of the input design and/or resulting gate sets can vary. In certain example embodiments, for instance, the process of FIG. 18 (or any of the processes disclosed herein) is used as part of a synthesis process generating a gate set description specifying 500 or more, 1000 or more, or 2500 or more gates. The particular operations or sequence should not be construed as limiting, as they can be performed alone or in any combination or subcombination with one another. Additionally, the illustrated operations can be performed together with one or more other operations. Still further, the identified operations need not be performed by a single software module, but can be implemented using multiple modules or software tools, which collectively perform the illustrated method.

At 1810, a target unitary described by a target angle and target precision is received (e.g., input, buffered into memory, or otherwise prepared for further processing).

At 1812, a corresponding quaternion approximation of the target unitary is determined. In particular embodiments, the corresponding quaternion approximation describes the target unitary in quaternion algebra. In some embodiments, and as discussed in detail above, the synthesis procedure can further comprise determining a cost vector for the target unitary (e.g., a cost vector describing a limit on a size of the single qubit gate set resulting from the synthesis procedure). In such embodiments, the act of determining the corresponding quaternion approximation can comprise finding one or more corresponding quaternion approximations that satisfy the cost vector. In some embodiments, the determining the corresponding quaternion approximation of the target unitary comprises selecting a first integer for use as a first quaternion element in the quaternion approximation, the first integer being selected so that a distance threshold to a R_(x), R_(y), or R_(z) rotation is satisfied. In such embodiments, a second integer for use as a second quaternion element in the quaternion approximation can also be selected, the second integer being selected so that, in combination with the first integer, the quaternion approximation produces a single qubit gate set that satisfies a cost threshold or gate set size limit. In particular implementations, the first integer and the second integer are selected from a ring of integers of the CM-field K. In particular implementations, the determining a corresponding quaternion approximation of the target unitary comprises finding a quaternion from the generalized Lipschitz order that has the following two properties: (1) d(U_(q),R_(z)(φ))≤ε; and (2) nrd(q)

_(F)=

₁ ^(L) ¹ . . .

_(M) ^(L) ^(M) , where d=√{square root over (1−|Tr(U_(q) ^(d)R_(z)(φ)|/2)},

₁, . . . ,

_(M) are the appropriate prime ideals of

_(F), and L₁, . . . , L_(M)∈

are their respective multiplicities in the decomposition of nrd(q)

_(F).

At 1814, the corresponding quaternion approximation is synthesized to produce a quantum circuit, the quantum circuit being over a single qubit gate set, get unitary, the single qubit gate set being realizable by a target quantum computer architecture. The single qubit gate set can be a Clifford+T basis gate set, a Clifford+e^(iπZ/12) basis gate set, a Clifford+e^(iπZ/16) basis gate set, a V-basis gate set, or any other single qubit base set expressible in terms of totally definite quaternion algebra.

At 1816, the quantum circuit is implemented in a target quantum circuit architecture. For instance, such implementation can be performed by a quantum circuit controller coupled to the target quantum circuit architecture and configured to implement the single qubit gate set for the target unitary in the target quantum circuit architecture. For example, embodiments as described in FIG. 21 can be used.

FIG. 19 is a flow chart 1900 showing another example synthesis process for producing a gate set in accordance with embodiments of the disclosed technology. The process can be performed, for instance, by a quantum computer synthesis tool adapted for use in a quantum computer design and implementation process and implemented by one or more computing devices. The size and complexity of the input design and/or resulting gate sets can vary. In certain example embodiments, for instance, the process of FIG. 19 (or any of the processes disclosed herein) generates a gate set description specifying 500 or more, 1000 or more, or 2500 or more gates. The particular operations or sequence should not be construed as limiting, as they can be performed alone or in any combination or subcombination with one another. Additionally, the illustrated operations can be performed together with one or more other operations. Still further, the identified operations need not be performed by a single software module, but can be implemented using multiple modules or software tools, which collectively perform the illustrated method.

At 1910, a program describing a desired computation to be performed in a target quantum computer architecture is input (e.g., buffered into memory or otherwise prepared for further processing).

At 1912, a gate set adapted for implementation on the target quantum computer architecture is generated based on the program. In the illustrated embodiment, the generating of the gate set includes generating a single qubit circuit for one or more respective unitaries used to perform the desired computation. Further, and as discussed in detail above, the generating the single qubit circuit comprises identifying a target unitary, and determining a corresponding quaternion approximation of the target unitary using a process in which one or more of the elements in the corresponding quaternion approximation are randomly selected.

At 1914, the gate set is implemented in a target quantum circuit architecture. For instance, such implementation can be performed by a quantum circuit controller coupled to the target quantum circuit architecture and configured to implement the gate set in the target quantum circuit architecture. For example, embodiments as described in FIG. 21 can be used.

In particular implementations, the determining the corresponding quaternion approximation of the target unitary comprises randomly selecting one or more values for use in the corresponding quaternion approximation from a constrained body of values that guarantees that a norm equation is solvable, and solving the norm equation to determine at least another one of the elements in the corresponding quaternion approximation. Further, in some implementations, the at least another one of the elements in the corresponding quaternion approximation is z, where z∈K, and wherein the norm equation for z is z(z*)=e, where e is an element of the totally real subfield F.

In some implementations, the determining the corresponding quaternion approximation of the target unitary comprises randomly sampling points from subsets of a ring of integers, and using the randomly selected sampling points as quaternion elements in the quaternion approximation. For example, the ring of integers can be a ring of integers of the CM-field K.

In certain implementations, the determining the corresponding quaternion approximation of the target comprises: partitioning a convex body of complex numbers into convex subsets, the partitioning being performed such that each subset has a lattice point corresponding to

_(K); randomly selecting one of the subsets; selecting a lattice point corresponding to

_(K) from the randomly selected subset; and using the selected lattice point as an element in the quaternion approximation.

FIG. 20 is a flow chart 2000 showing another example synthesis process for producing a gate set in accordance with embodiments of the disclosed technology. The process can be performed, for instance, by a quantum computer synthesis tool adapted for use in a quantum computer design and implementation process and implemented by one or more computing devices. Although the process in FIG. 20 is described with respect to a single qubit gate sets, the size and complexity of the input design and/or resulting gate sets can vary. In certain example embodiments, for instance, the process of FIG. 20 (or any of the processes disclosed herein) is used as part of a synthesis process generating a gate set description specifying 500 or more, 1000 or more, or 2500 or more gates. The particular operations or sequence should not be construed as limiting, as they can be performed alone or in any combination or subcombination with one another. Additionally, the illustrated operations can be performed together with one or more other operations. Still further, the identified operations need not be performed by a single software module, but can be implemented using multiple modules or software tools, which collectively perform the illustrated method.

At 2010, a target unitary described by a target angle and target precision is received (e.g., input, buffered into memory, or other prepared for further processing).

At 2012, a corresponding quaternion approximation of the target unitary is determined.

At 2014, the corresponding quaternion approximation is synthesized to produce a single qubit circuit for the target unitary, the single qubit circuit being realizable by a quantum computer architecture.

Further, in the illustrated embodiment and as discussed in detail above, the act of determining the corresponding quaternion approximation of the target unitary uses a Closest Vector Problem (CVP) technique to select values for the corresponding quaternion approximation. In some embodiments, the determining the corresponding quaternion approximation of the target unitary comprises selecting values from a Hermite-Korkine-Zolotarev, Block-Korkine-Zolotarev, Lenstra-Lenstra-Lovasz, or other size-reduced basis as an element in the corresponding quaternion approximation. In certain embodiments, the determining the corresponding quaternion approximation of the target unitary comprises using a Lenstra-Lenstra-Lovasz or other lattice basis reduction technique.

At 2016, the single qubit circuit for the target unitary is implemented in a target quantum circuit architecture. For instance, such implementation can be performed by a quantum circuit controller coupled to the target quantum circuit architecture and configured to implement the single qubit gate set for the target unitary in the target quantum circuit architecture. For example, embodiments as described in FIG. 21 can be used.

10. Representative Computing Environments

FIG. 21 and the following discussion are intended to provide a brief, general description of an exemplary computing environment in which the disclosed technology may be implemented. Although not required, the disclosed technology is described in the general context of computer executable instructions, such as program modules, being executed by a personal computer (PC). Generally, program modules include routines, programs, objects, components, data structures, etc., that perform particular tasks or implement particular abstract data types. Moreover, the disclosed technology may be implemented with other computer system configurations, including hand held devices, multiprocessor systems, microprocessor-based or programmable consumer electronics, network PCs, minicomputers, mainframe computers, and the like. The disclosed technology may also be practiced in distributed computing environments where tasks are performed by remote processing devices that are linked through a communications network. In a distributed computing environment, program modules may be located in both local and remote memory storage devices.

With reference to FIG. 21, an exemplary system for implementing the disclosed technology includes a general purpose computing device in the form of an exemplary conventional PC 2100, including one or more processing units 2102, a system memory 2104, and a system bus 2106 that couples various system components including the system memory 2104 to the one or more processing units 2102. The system bus 2106 may be any of several types of bus structures including a memory bus or memory controller, a peripheral bus, and/or a local bus using any of a variety of bus architectures. The exemplary system memory 2104 includes read only memory (ROM) 2108 and random access memory (RAM) 2110. A basic input/output system (BIOS) 2112, containing the basic routines that help with the transfer of information between elements within the PC 2100, is stored in ROM 2108. As shown in FIG. 21, the system memory 2104 stores computer-executable instructions for performing any of the disclosed techniques (e.g., determining quaternion approximations, quaternion mappings, defining special unitaries, finding gate sets, algebraic computations, and exact/approximate circuit synthesis) in respective memory portions (shown generally as executable software 2116 for performing any embodiment of the disclosed synthesis techniques)).

The exemplary PC 2100 further includes one or more storage devices 2130, such as a hard disk drive for reading from and writing to a hard disk, a magnetic disk drive for reading from or writing to a removable magnetic disk, and/OR an optical disk drive for reading from or writing to a removable optical disk (such as a CD-ROM or other optical media). Such storage devices can be connected to the system bus 2106 by a hard disk drive interface, a magnetic disk drive interface, and/or an optical drive interface, respectively. The drives and their associated computer readable media provide nonvolatile storage of computer-readable instructions, data structures, program modules, and other data for the PC 2100. Other types of computer-readable media which can store data that is accessible by a PC, such as magnetic cassettes, flash memory, digital video disks, CDs, DVDs, RAMs, NVRAMs, ROMs, and the like, may also be used in the exemplary operating environment. As used herein, the terms storage, memory, and computer-readable media do not include or encompass propagating carrier waves or signals per se.

A number of program modules may be stored in the storage devices 2130, including an operating system, one or more application programs, other program modules, and program data. Storage of results of quantum syntheses and instructions for obtaining such syntheses (e.g., instructions for performing any embodiment of the disclosed technology) can be stored in the storage devices 2130. A user may enter commands and information into the PC 2100 through one or more input devices 2140 such as a keyboard and a pointing device such as a mouse. Other input devices may include a digital camera, microphone, joystick, game pad, satellite dish, scanner, or the like. These and other input devices are often connected to the one or more processing units 2102 through a serial port interface that is coupled to the system bus 2106, but may be connected by other interfaces such as a parallel port, game port, or universal serial bus (USB). A monitor 2146 or other type of display device is also connected to the system bus 2106 via an interface, such as a video adapter. Other peripheral output devices, such as speakers and printers (not shown), may be included. In some cases, a user interface is displayed so that a user can input a circuit for synthesis, and verify successful synthesis.

The PC 2100 may operate in a networked environment using logical connections to one or more remote computers, such as a remote computer 2160. In some examples, one or more network or communication connections 2150 are included. The remote computer 2160 may be another PC, a server, a router, a network PC, or a peer device or other common network node, and typically includes many or all of the elements described above relative to the PC 2100, although only a memory storage device 2162 has been illustrated in FIG. 21. The personal computer 2100 and/or the remote computer 2160 can be connected to a logical a local area network (LAN) and a wide area network (WAN). Such networking environments are commonplace in offices, enterprise wide computer networks, intranets, and the Internet.

When used in a LAN networking environment, the PC 2100 is connected to the LAN through a network interface. When used in a WAN networking environment, the PC 2100 typically includes a modem or other means for establishing communications over the WAN, such as the Internet. In a networked environment, program modules depicted relative to the personal computer 2100, or portions thereof, may be stored in the remote memory storage device or other locations on the LAN or WAN. The network connections shown are exemplary, and other means of establishing a communications link between the computers may be used.

With reference to FIG. 22, an exemplary system for implementing the disclosed technology includes computing environment 2200, where synthesis (compilation) into implementable gate sets (e.g., using any of the disclosed embodiments) is separated from the quantum processing that consumes the compiled circuits. The environment includes one or more quantum processing unit(s) 2202 and one or more monitoring/measuring device(s) 2246. The quantum processing unit(s) execute quantum circuits that are synthesized by synthesis/compiler unit 2220 (e.g., configured to perform any of the disclosed procedures) using one or more classical processor(s) 2210. The synthesized/compiled quantum circuits are downloaded into or used to program or configure the quantum processing unit(s) 2202 (e.g., via control lines (quantum bus) 2206). Synthesis procedures according to any of the disclosed embodiments (e.g., synthesis procedures associated with quaternion approximations and gate sets are stored in a memory 2221, and special unitaries or procedures associated with such unitaries are stored in a memory 2203.

With reference to FIG. 22, synthesis (or compilation) is the process of translating a high-level description of a quantum algorithm into gate sets (e.g., a sequence of quantum circuits. Such high-level descriptions may be stored, as the case may be, on one or more external computer(s) 2260 outside the computing environment 2200 utilizing one or more memory and/or storage device(s) 2262, then downloaded as necessary into the computing environment 2200 via one or more communication connection(s) 2250. Gate sets (circuits) 2206 defined using gate sets from quaternion algebra based synthesis (according to any of the disclosed embodiments) are coupled to the quantum processor 2202.

The quantum processing unit(s) can be one or more of, but are not limited to: (a) a superconducting quantum computer; (b) an ion trap quantum computer; or (c) a fault-tolerant architecture for quantum computing (e.g., a topological quantum computer using Majorana zero modes). The synthesized gate sets (e.g., using any of the disclosed embodiments) can be sent into (or otherwise applied to) the quantum processing unit(s) via control lines 2206 at the control of the classical processor 2210 and/or synthesis unit 2220. In the illustrated example, the desired quantum computing process is implemented with the aid of one or more QP subcontrollers 2205 that are specially adapted to control a corresponding one of the quantum processor(s) 2202. For instance, in one example, the classical processor 2210 and/or synthesis unit 2220 facilitates implementation of the compiled quantum circuit by sending instructions to one or more memories (e.g., lower-temperature memories), which then pass the instructions to low-temperature control unit(s) (e.g., QP subcontroller(s) 2205) that transmit, for instance, pulse sequences representing the gates to the quantum processing unit(s) 2202 for implementation. In other examples, the QP subcontroller(s) 2205 operate to provide appropriate magnetic fields, encoded operations, or other such control signals to the quantum processing unit(s) 2202 to implement the operations of the compiled quantum computer circuit description. The classical processor 2210 can further interact with measuring/monitoring devices (e.g., readout devices) 2246 to help control and implement the desired quantum computing process (e.g., by reading or measuring out data results from the quantum processing units once available, etc.)

Having described and illustrated the principles of the disclosed technology with reference to the illustrated embodiments, it will be recognized that the illustrated embodiments can be modified in arrangement and detail without departing from such principles. For instance, elements of the illustrated embodiments shown in software may be implemented in hardware and vice-versa. Also, the technologies from any example can be combined with the technologies described in any one or more of the other examples. It will be appreciated that procedures and functions such as those described with reference to the illustrated examples can be implemented in a single hardware or software module, or separate modules can be provided. The particular arrangements above are provided for convenient illustration, and other arrangements can be used. 

What is claimed is:
 1. A quantum circuit synthesizer system, comprising: a processor; and at least one memory coupled to the processor and having stored thereon processor-executable instructions for a quantum computer synthesis procedure that comprises: receiving a target unitary described by a target angle and target precision; determining a corresponding quaternion approximation of the target unitary; and synthesizing a quantum circuit corresponding to the quaternion approximation, the circuit being over a single qubit gate set, the single qubit gate set being realizable by a target quantum computer architecture, wherein the determining a corresponding quaternion approximation of the target unitary comprises finding a quaternion from an order of a totally definite quaternion algebra defined over a totally real number field F that has the following two properties: d(U _(q) ,R _(z)(φ))≤ε;  (1) and nrd(q)

_(F)=

₁ ^(L) ¹ . . .

_(M) ^(L) ^(M) ,  (2) wherein d is a distance function,

₁, . . . ,

_(M) are the appropriate prime ideals of

_(F) ₁ , L₁, . . . , L_(M)∈

are their respective multiplicities in the decomposition of nrd(q)

_(F), and nrd(q) is a reduced norm of a quaternion q.
 2. The quantum circuit synthesizer system of claim 1, wherein the synthesis procedure further comprises determining a cost vector for the target unitary, and wherein the determining the corresponding quaternion approximation comprises finding one or more corresponding quaternion approximations that satisfy the cost vector, wherein the cost vector describes a limit on a size of the single qubit gate set resulting from the synthesis procedure.
 3. The quantum circuit synthesizer system of claim 1, wherein the distance function d is d=√{square root over (1−|Tr(U _(q) ^(d) R _(z)(φ))|/2)}.
 4. The quantum circuit synthesizer system of claim 1, wherein the single qubit gate set is a Clifford+e^(iπZ/12) basis gate set, a Clifford+e^(iπZ/16) basis gate set, a V-basis gate set, or any other single qubit base set expressible in terms of totally definite quaternion algebra.
 5. The quantum circuit synthesizer system of claim 1, further comprising a quantum circuit controller coupled to the target quantum circuit architecture and configured to implement the quantum circuit in the target quantum circuit architecture.
 6. A quantum circuit synthesizer system, comprising: a processor; and at least one memory coupled to the processor and having stored thereon processor-executable instructions for a quantum computer synthesis procedure that comprises: receiving a target unitary described by a target angle and target precision; determining a corresponding quaternion approximation of the target unitary; and synthesizing a quantum circuit corresponding to the quaternion approximation, the circuit being over a single qubit gate set, the single qubit gate set being realizable by a target quantum computer architecture, wherein the determining the corresponding quaternion approximation of the target unitary comprises: selecting a first algebraic integer for use as a first quaternion element in the quaternion approximation, the first algebraic integer being selected so that a distance threshold to a Rx, Ry, or Rz rotation is satisfied.
 7. The quantum circuit synthesizer system of claim 6, wherein the determining the corresponding quaternion approximation of the target unitary further comprises: selecting a second algebraic integer for use as a second quaternion element in the quaternion approximation, the second algebraic integer being selected so that, in combination with the first algebraic integer, the quaternion approximation produces a quaternion that satisfies a cost limit and the corresponding unitary satisfies the distance threshold.
 8. A method, comprising: by a quantum computer synthesis tool adapted for use in a quantum computer design and implementation process and implemented by one or more computing devices: inputting a program describing a desired computation to be performed in a target quantum computer architecture; generating a gate set adapted for implementation on the target quantum  computer architecture based on the program, wherein the generating includes generating a single qubit circuit for one or more respective unitaries used to perform the desired computation, wherein the generating the single qubit circuit comprises: identifying a target unitary, and determining a corresponding quaternion approximation of the target unitary using a process in which one or more of the elements in the corresponding quaternion approximation are randomly selected, wherein the determining the corresponding quaternion approximation of the target unitary comprises: randomly selecting one or more values for use in the corresponding quaternion approximation from a constrained body of values that guarantees that a norm equation is solvable; and solving the norm equation to determine at least another one of the elements in the corresponding quaternion approximation.
 9. The method of claim 8, wherein the at least another one of the elements in the corresponding quaternion approximation is z, where z is an algebraic integer from K and wherein the norm equation for Z is Z(Z*)=e₁, where e is an element of the totally real subfield F.
 10. A method, comprising: by a quantum computer synthesis tool adapted for use in a quantum computer design and implementation process and implemented by one or more computing devices: inputting a program describing a desired computation to be performed in a target quantum computer architecture; generating a gate set adapted for implementation on the target quantum  computer architecture based on the program, wherein the generating includes generating a single qubit circuit for one or more respective unitaries used to perform the desired computation, wherein the generating the single qubit circuit comprises: identifying a target unitary, and determining a corresponding quaternion approximation of the target unitary using a process in which one or more of the elements in the corresponding quaternion approximation are randomly selected, wherein the determining the corresponding quaternion approximation of the target unitary comprises: randomly sampling points from subsets of a ring of integers; and using the randomly selected sampling points as quaternion elements in the quaternion approximation.
 11. A method, comprising: by a quantum computer synthesis tool adapted for use in a quantum computer design and implementation process and implemented by one or more computing devices: inputting a program describing a desired computation to be performed in a target quantum computer architecture; generating a gate set adapted for implementation on the target quantum  computer architecture based on the program, wherein the generating includes generating a single qubit circuit for one or more respective unitaries used to perform the desired computation, wherein the generating the single qubit circuit comprises: identifying a target unitary, and determining a corresponding quaternion approximation of the target unitary using a process in which one or more of the elements in the corresponding quaternion approximation are randomly selected, wherein the determining the corresponding quaternion approximation of the target comprises: partitioning a convex body of complex numbers into convex subsets, the partitioning being performed such that each subset has a lattice point corresponding to

_(K); randomly selecting one of the subsets; selecting a lattice point corresponding to

_(K) from the randomly selected subset; and using the selected lattice point as an element in the quaternion approximation.
 12. One or more computer-readable memory or storage devices storing computer-executable instructions which when executed by a computer cause the computer to perform a quantum computer synthesis procedure comprising: receiving a target unitary described by a target angle and target precision; determining a corresponding quaternion approximation of the target unitary; synthesizing the corresponding quaternion approximation to produce a single qubit circuit for the target unitary, the single qubit circuit being realizable by a quantum computer architecture, wherein the determining the corresponding quaternion approximation of the target unitary uses a Closest Vector Problem (CVP) technique to select values for the corresponding quaternion approximation.
 13. The one or more computer-readable memory or storage devices of claim 12, wherein the determining the corresponding quaternion approximation of the target unitary comprises selecting values from a size-reduced basis as an element in the corresponding quaternion approximation. 